「成功するまでの回数」の期待値を求めるCプログラムとRubyスクリプト

いきなりですが問題です．

成功確率が $p$ (ただし[tex:0

直感でこの期待値は $1/p$ ，すなわち逆数になるのではないかと考えます． $p=1/2$ なら，期待値は2，すなわち平均的に見れば，2回で終了する，ということです． $p=1/10$ なら，期待値は10です．もちろん1回で終了するケースも，何十回やっても終わってくれないケースもあるでしょうが，「平均的に」見たときの回数を教えてくれるのが，期待値というものです．
この期待値を数式で表すと，
$1\times p + 2\times(p-1)p + \cdots + n\times (p-1)^{n-1}p + \cdots$
すなわち
$\sum_{n=1}^{\infty}n(p-1)^{n-1}p$
になります．ここから式変形に奮闘しましたが，どうも思っていた解になってくれません．計算ミスだと思うthere must be an error in calculationのですが，どこか分からない…．
ここでシミュレーションに頼ってみました．まずはRubyで，それからCで書きました*1．

/*
 * expectation.c
 *   独立試行で「1回成功すればそこで終了」の回数の期待値を
 *   シミュレーションにより求める．
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/* 期待値を求める */
void calc_expectation(int pinv, int exec)
{
  int i;
  double sum;  /* 試行回数の和 */

  /* 引数チェック(0除算をしないように) */
  if (pinv == 0) {
    printf("p = infinity?\n");
    return;
  }
  if (exec == 0) {
    printf("execution = 0\n");
    return;
  }

  for (i = 1; i <= exec; i++) {
    double count = 1; /* 成功するまでの試行回数 */

    for (count = 1; (unsigned long int)random() % pinv != 0; count++) {
      printf(".");
    }

    printf("! %d\n", (int)count);

    sum += count;
  }

  printf("\n");
  printf("p = 1/%d\n", pinv);
  printf("execution = %d\n", exec);
  printf("expectation = %f\n", sum / exec);
}

int main(int argc, char *argv[])
{
  int pinv = 2; /* 独立試行の成功確率の逆数 */
  int exec = 10000; /* 実行回数 */

  if (argc >= 2) {
    pinv = atoi(argv[1]);
  }
  if (argc >= 3) {
    exec = atoi(argv[2]);
  }

  srandom(time(NULL));

  calc_expectation(pinv, exec);
  
  return 0;
}

#!/usr/bin/env ruby

# expectation.rb
#   独立試行で「1回成功すればそこで終了」の回数の期待値を
#   シミュレーションにより求める．

class Expectation
  def initialize(pinv_ = 2, exec_ = 10000)
    @pinv = pinv_.to_i  # 独立試行の成功確率の逆数
    @exec = exec_.to_i  # 実行回数
    @verbose = true  # シミュレーション状況を出力するならtrue
  end

  def try
    # 1回成功するまで独立試行を行い，その回数を返す．

    count = 0
    while true
      count += 1
      if rand(@pinv) == 0
        puts "! #{count}"  if @verbose
        return count
      end
      print "."  if @verbose
    end
  end

  def start
    sum = 0  # 試行回数の和
    1.upto(@exec) do
      sum += try
    end

    e = sum.to_f / @exec.to_f  # 期待値

    puts ""
    puts "p = 1/#{@pinv}"
    puts "execution = #{@exec}"
    puts "expectation = #{e}"
  end
end

if __FILE__ == $0
  pinv = ARGV.shift || 2  # 独立試行の成功確率の逆数
  exec = ARGV.shift || 10000  # 実行回数
  Expectation.new(pinv, exec).start
end

とることのできるコマンドライン引数は，どちらも同じです．最初の引数は，成功確率の逆数で，整数で与えるものとします．2番目の引数は，「成功するまでその試行を繰り返」すのを1回の実行として，何回実行するかで，これももちろ整数です．コマンドライン引数を書かずに実行したときの値，すなわちデフォルト値は，それぞれ2と10000です．
冒頭で挙げたとおり，ベルヌーイ試行での生起確率は0より大きく1未満の値ですが，ここでの入力ではその逆数をとり，プログラムではpinvという変数(「p inverse」の略です)に格納します．
「成功判定」は，CとRubyとで若干違います．Cでは，randomの戻り値はlong int型の整数値なので，pinvで割ってその余りが0なら*2，成功としました．負の数を割らないよう，unsignedつきのキャストもしています．一方Rubyでは，「rand(@pinv)」と書けば，0以上@pinv未満の整数を返してくれるので，単純にその値を0と比較し，等しいなら成功としました．
さて実行．期待値はpinvになってくれるといいのですが，回数が多くして実行すると，確かにpinv，すなわち生起確率の逆数に近づいています．
…と，プログラミングに時間を費やしてから，計算していたメモに目をやると，間違いに気付き，無事に $1/p$ という結果に至りました．
ところでなぜこんな問題を考えたのかというと，以前に

1024ビットの乱数で，1のビットがちょうど512個になる確率は，40分の1程度．思っていたより高いなあというのが，率直なところです．
半分ビットの確率を求めるRubyスクリプト - わさっき

と書いたのですが，1024ビットの乱数を生成して，1のビットがちょうど512個になるまで生成するとしたとき，生成回数の期待値を確認したかったからでした．

*1:Rubyでは，initializeメソッドの変数@verboseにtrueを書いていますが，ここをfalseにすることによって，途中経過の出力を抑制します．Cでは面倒なので書いていません．途中経過は，試行に失敗したら「.」，成功したら「!」と成功までの試行回数を出力するというものです．

*2:最初に「unsigned long int型の最大値」をpinvで割った値を閾値として新たに変数に格納しておき，randomで得られる値をこれと比較するようにすれば，乗算をしない分，速度が速くなりますね．

わさっきhb

大学（教育研究）とか，親馬鹿とか，和歌山とか，とか，とか．

「成功するまでの回数」の期待値を求めるCプログラムとRubyスクリプト