統計・AI

manao55.hatenablog.com

manao55.hatenablog.com

正規乱数表現

正規乱数

ja.wikipedia.org

を平均, を標準偏差とすると, 正規分布 は次のように表現される分布となります。

画像引用元：「Metrology: The Science of Measurement」https://www.muelaner.com/wp-content/uploads/2013/07/Standard_deviation_diagram.png

この平均 , 標準偏差 の正規分布 に従う乱数は, ボックス・ミュラー法により, 2個の一様乱数 , を用いて

で得ることができます。

プログラム

#include <iostream>  /* cout */
#include <math.h>    /* sin, cos */
#include <cstdlib>   /* rand */
#include <iomanip>   /* setw */
using namespace std;

void boxrnd(double, double, double *, double *);

int main() {
    int i, j, hist[1000];
    double x,y;

    for (i=0;i<1000;i++)
        hist[i] = 0;

    for (i=0;i<10000;i++) {
        boxrnd(2.6,10.0,&x,&y);
        hist[(int)x]++;
        hist[(int)y]++;
    }

    for (i=0;i<20;i++){
        cout << setw(3) << i+1 << " : " << setw(4) << hist[i] << " I " ;
        for(j=1;j<=hist[i]/100;j++){
            cout << "*" ;
        }
        cout << endl;
    }
    return 0;
}

void boxrnd(double sig, double m, double *x, double *y){

    double r1,r2;
    r1 = rand()/32767.1;
    r2 = rand()/32767.1;
    *x = sig*sqrt(-2*log(r1))*cos(2*3.14159*r2) + m;
    *y = sig*sqrt(-2*log(r1))*sin(2*3.14159*r2) + m;
}

実行結果

  1 :    5 I
  2 :   17 I
  3 :   50 I
  4 :  141 I *
  5 :  335 I ***
  6 :  730 I *******
  7 : 1267 I ************
  8 : 1939 I *******************
  9 : 2613 I **************************
 10 : 3069 I ******************************
 11 : 2987 I *****************************
 12 : 2528 I *************************
 13 : 1897 I ******************
 14 : 1207 I ************
 15 :  669 I ******
 16 :  333 I ***
 17 :  140 I *
 18 :   47 I
 19 :   19 I
 20 :    6 I

繰り返し制御の紹介

はじめに

最近, 機械学習, 深層学習の発展や社会実装に伴い, 数多くのAIベンチャーが目に移ります。主な適用先は, 自動運転車や産業用ロボット, バイオ分野で, 特に画像認識の応用が見受けられます。

ロボット・AIというワードをみると, 制御と学習の関連が連想されます。 ひとによって異なるのでしょうが, AIで動くロボットという概念はどこかにいる人類が考えていそうな欲のひとつなのかなと思われます。（ロボットが自律して判断して行動すると面白そうですよね）

制御と学習は昔から思いつく人は多いでしょうし, 始まりはなんだったのだろうかとふと思いました。近頃のロボティクスなどの研究では比較的新しい学習理論が適用されるとは思いますが, ここでは温故知新ということで, 個人的に学習と制御に関する初期の理論を調べてみて, 繰り返し制御のアイデアが興味を引いたので取り上げて確認しようと思いました。

繰り返し制御

繰り返し制御は, 東京工業大学制御工学科の井上 悳氏, 中野 道雄氏 およびに高エネルギー研 (現高エネルギー加速器研究機構）の久保 忠志氏, 松本 啓氏, 馬場 斉氏らにより陽子シンクロトロン加速器の電源の制御のために最初に考案されました。

井上 悳, 中野 道雄, 久保 忠志, 松本 啓, 馬場 斉, 陽子シンクロトロン電磁石電源の繰返し運転における高精度制御, 電気学会論文誌C, 1980, 100 巻, 7 号, p. 234-240 https://www.jstage.jst.go.jp/article/ieejeiss1972/100/7/100_7_234/_article/-char/ja

1周期 (試行) 前の偏差を利用して周期的な目標入力に高精度で追従させる制御系のようです。 当時は理論的に, またさまざまな分野への応用に関しての研究がなされていたようです。

原 辰次, 小俣 透, 中野 道雄, 繰返し制御系の安定条件と設計法, 計測自動制御学会論文集, 22 巻, 1 号, p.36-42 https://www.jstage.jst.go.jp/article/sicetr1965/22/1/22_1_36/_pdf

小林 史典, 原 辰次, 田中 啓友, 中野 道雄, 繰返し制御を応用したモータの回転むら低減法, 電気学会論文誌D, 107 巻, 1 号, p.29-34 https://www.jstage.jst.go.jp/article/ieejias1987/107/1/107_1_29/_pdf

ずいぶん昔の理論ではありますが, 理論的な考察も多く, 実応用に期待されていた制御法であったことは伺えます。採用するメリットとして,

・ロボットの繰返し動作など, 周期的な目標入力に対して有効である. （閉ループ内に周期信号の発生機構をもつ繰返し制御系は,電源周波数やモータの回転 数に依存した周期的な外乱入力の除去にも効果がある.）

・逆システムの反復的生成法と考えられる学習制御方式は,ロボットの軌道制御などを実現する制御入力の自動生成の強力なひとつの方式である.

・数少ない日本人考案の国産の制御法である.

大元は加速器がターゲットでしたが, 1970年代ごろ複数の日本の企業が産業用ロボットの生産を開始したことにより, 反復作業の高性能化を狙った制御手法を求めていて, 学習制御のひとつとして注目され始めたのかなと推測しています。

加速器を安定に制御させて得たかった彼らの本当の目的についてはここでは触れないことにします。

繰り返し制御の構成

繰り返し制御の構成方法として様々な方法がありますが, 基本的な制御ブロックは以下のようになります。

または のブロックをループ内の前向き伝達関数のブロックに配置して以下のようにもなります。

制御器 の前段に設置されることが多いです。 二自由度制御のフィードフォワード制御部として扱われることがあります。

原理

繰り返し制御の原理について説明します。定式化の都合上, いくつか前提条件を設定します。

(ⅰ) 制御系の目標値 は一定周期 を持つ繰返し信号か, あるいは一定値信号であるとする｡

(ⅱ) この目標値のもとで, 波高値の有限な何らかの操作信号を制御対象に加えることにより制御偏差をなくすることが可能であるとする。

(ⅲ) 制御対象は線形であるとする｡

前章のブロック図で示すような繰り返し制御器の伝達関数を とすると,

周波数応答は, として,

の仮定条件およびに上式より, は離散角周波数

において各々一定の周波数成分を持ちます。（極を持ちます）

このとき, 周波数成分は

となり, 利得が無限大に近づくことが分かります。

以下に繰り返し制御部の周波数応答を示します。

% repete_bode.m
s = tf('s');
T = 2;    % 2秒の周期関数の発生機構
G = exp(-s*T)/(1-exp(-s*T));
bode(G);

このため, 制御系が安定で, が有限の周波数で零にならず, 外乱 が一定値あるいは の整数分の１の周期を持つとき、偏差は零へ収束します。 例えば外乱がないとき, 閉ループ伝達関数 は,

1 に収束するので, 偏差が零へと収束していきます。

サーボ系の重要な要件である内部モデル原理 (前向き伝達関数のいずれかに目標値の伝達関数が含まれれば定常偏差が0になる) がこの場合でも成立するのであれば, 周期関数の発生機構を制御ループ内に挿入することで, 周期 の任意の周期目標入力に対して定常偏差をなくすことができます。

強引な式変形ですが,

と近似すると,

と表してみると, 周期関数の発生機構自体が（ゆっくりとした）近似的な積分器として解釈することもできます。

コントローラ やコントローラ などをかんがみると, 1形の制御系と繰り返し制御系とは

といった対応があることが分かります。

以下の図において, 単純に周期関数発生機構 だけでなく を並列に接続し, PI制御器を模したような繰り返し補償器としての構成も考えることができます。この構成を, 修正繰り返し制御系とも呼びます。

初めの繰り返し制御のブロック図は としたときに対応することとなります。

PI制御器と繰り返し制御器はよく似ていますが, 前者が有限次元システムであるのに対して, 繰り返し制御はむだ時間 をもつ無限次元システムである点が異なります。（次数を問いません） この性質は制御の安定性を議論する際に特殊な側面が表れてきます。

繰り返し制御の安定条件

初めの図の制御系におきまして

が成立します。（, はそれぞれ, の初期値応答のラプラス変換である。） これらの関係を用いると, , 式より

が成立します。ここで, , , 式より,

が導かれ, にを代入すると,

つまり,

となります。ここで,

となります。

式の関係をブロック線図で表すと, 以下の等価系が得られます。

ではこの図の等価系の入出力条件から、繰り返し制御系の安定条件を求めてみましょう。

とおいて

と定義します。 で , が有界であるならば, , となります。

は

でノルムが定義される ]上の加速関数 の集合です。

ここで, ^{-1} ] の漸近安定性（ で平衡点 に収束）を仮定すると, 式より,

となります。等価系の閉ループにおいてスモールゲイン定理を用いると,

が成立します。

, の仮定の下, 有界で連続な周期信号 , に対して

が成立し, 繰り返し制御の定理となります。

式の考察について, 式は,

と変形することができ,

と置くと, 以下のように場合分けすることができます。

> ( > )のとき,

= ( = ) のとき,

< ( < )のとき,

と導くことができます。特に のとき, 最適制御やカルマンフィルタの最適性条件 (円条件) と等しくなります。以下にナイキスト線図を示します.（すみません.手書きですが）

繰り返し制御の設計法

当時, 繰り返し制御などの学習制御の収束条件についても最適制御の最適性条件との関係性が考察されており, 状態フィードバックによる設計法が提案されています。

川村 貞夫, 宮崎 文夫, 有本 卓, 学習制御システムのシステム論的考察, 計測自動制御学会論文集, 21 巻, 5 号, p.445-450 https://www.jstage.jst.go.jp/article/sicetr1965/21/5/21_5_445/_pdf/-char/ja

個人的にはあまり本質的な場面ではない（は？）と思ってしまったので, ここでは出力動的フィードバックによる繰返し制御系の設計法をなるべく簡単に事実だけ述べて説明します。

下図の制御系において,

プラント が与えられたとき, が

の二つの条件を満たす補償器 , の設計法について検討します。

ここではカルマンフィルタと完全制御の手法を用いた最小位相系に対する修正繰り返し制御の設計法 を説明します。

まず, を最小実現 を用いて,

と置きます。このとき, を 以下の図のように構成します。

このとき, は,

と表せます.　ここで, を完全制御のアルゴリズム†（補足）で求め, それを を表すと]

が成立します。これを 式に適用すると,

このとき, はカルマンフィルタのゲイン

: 非負定, : 可観測対という条件で, はリカッチ方程式

で定めるとすると, ロバスト性不等式

つまり

が成立します。したがって,

と設定すれば, で は を除いて漸近的に式 を満たします。 この時, 式より,

であるので, カルマンフィルタのゲイン減少に関する安定余裕が であることを用いると, 式 式が成立することが確かめられます。 これまでの議論をまとめると, 漸近的設計法は以下のようにまとめられます。

の最小実現 を求める. を図のように構成し, ゲイン と, を以下のように決める。

: カルマンフィルタのゲインの 倍

: 完全制御のゲイン

式 を満たすように , を決定する。

あとがき

正直自分で書いていてこれだけでは自分がよくわかってないなと思ったので(???), 具体的なプラントを使って制御する様子を Matlab/Simulink advent Calender 4日目の記事で書いてみました。ご興味ありましたらご覧ください。

外乱オブザーバと繰り返し制御を併用した単軸アームの学習制御 https://qiita.com/Manao/items/cb8c0ab05cee98e6fde3

（補足）完全制御について

線形システム

を考える。ここで は 次元状態ベクトル, は, 次元入力ベクトル, は 次元出力ベクトルで, , , は適当なサイズの定数行列であります。上システムは、可制御・可観測であるとします。評価関数

に対する最適レギュレータは, 状態フィードバック

で与えられます。ただし　 はリカッチ方程式

の正定解です。最適還送差行列を

で定義すると, は不等式

を満足することが知られています。上式は安定余裕や感度改善などの最適レギュレータのロバスト性に関する議論の基礎を与えるもので, ロバスト性不等式と呼ばれます。

カルマンフィルタのロバスト性不等式も双対性を用いれば導かれます。つまり,

で与えられるカルマンフィルタゲインに対して, 最適還送差行列を

で定義すると, は

を満足します。ここで, , はそれぞれ外乱, 観測雑音の共分散行列である。 完全制御について, 必要な結果を述べます。パラメータ を含む状態フィードバックの族

を考えます。ただし はゲイン行列の各要素に有利関数として含まれる。 が

を満足するとき, はシステムに対して完全制御を達成するといいます。 以下 で閉ループ系の極が次の条件を満足する場合を考えます。

[A] のすべての固有値（状態フィードバックのもとでの閉ループ極）は次の２つの条件のうちいずれかを満足する。

ただし は負の実部を持つ複素数である。仮定[A]のもとで関係(20)は等価な周波数値域における条件

に置き換えられます。完全制御の存在条件は次のように与えられます。

・命題1

が存在するための必要十分条件は, システム が右可逆で最小位相系（すべての零点の実部が負）であることです。

完全制御の双対は完全観測です。完全観測は を満足するオブザーバゲインの族 によって定義される。つまり恒等オブザーバ

は, (20)式に対応した関係

を満足します。[A]と同様の仮定の下で (21)式は周波数領域における条件

と等価になります。完全観測を と書きます。完全観測の存在は[命題1]で右可逆を左か逆に置き換えたものであります。

外乱オブザーバと繰り返し制御を併用した単軸アームの学習制御

制御性能の向上を目指す

manao55.hatenablog.com

上の記事の続きとなります。

前回では, 制御系の設計, 外乱オブザーバの設計を行い, 外乱が発生してもアームが安定化する制御を行いました。

しかし定常的には安定する制御ではあったものの, 過渡の部分でオーバーシュートが発生していました。

今回は制御性能を上げることを目的に, 繰り返し制御を用いて過渡応答を向上させる制御を実装したいと思います。今回の"学習"は繰り返し制御の特長である, 「一定時間前の試行を利用して制御を向上させること」, と定義させていただきます。

では, よろしくお願いします。

制御ブロック

制御ブロックについて説明します。前回では外乱オブザーバと PD 制御器を併用した構成となっておりました。オレンジ色の部分がオブザーバとなっております。

今回は, 繰り返し制御を追加したため制御ブロックは以下のようになります。

緑色で囲った部分が追加した繰り返し制御部分となります。 この制御を加えることによって, 制御性能が向上すると望ましいのですが..

以下にSimulink モデルの図（デジタル制御）を示します。

繰り返し制御の実装

自分で試行錯誤で構成した繰り返し制御系を作ってみました。おそらく修正繰り返し制御だと思います。

今回はマイコンを使用したデジタル制御にて実装を行うので, 離散時間系列で考えるため,

といった変換を行います。 は周期（今回は4秒）であり, はサンプリング数の周期なので, 今回は,

となります。

偏差 の部分に一周期前の偏差 を追加するような制御で補正しつつ, 前向き伝達関数の比例ゲイン (学習係数) で補正量を調整できると解釈しています。

システム同定にてパラメータを取得しておきます。

以下に実装したプログラムを示します。

% pos_pd_mbd.m
s = tf('s');
N = 200;

% Initialize & load data
close all
clear all
load sim_param
load model_data

% PD design by pole placement
omega_n = 12;
zeta    = 0.8;
p1 = (-zeta + j*sqrt(1-zeta^2))*omega_n;
p2 = (-zeta - j*sqrt(1-zeta^2))*omega_n;

% Set PD parameters
Kp = p1*p2*T/K;
Kd = -((p1+p2)*T + 1)/K;
Ki = 0;

% Set learning gain
Kl = 0.15;

% Display PI parameters
disp('>>> PI parameters <<<')
fprintf('Kp  = %f\n',Kp);
fprintf('Ki  = %f\n',Ki);
fprintf('Kd  = %f\n',Kd);

% observer setting
omega_c = 7.5;        % cutoff frequency
Fs = omega_c^2 / (s^2 + sqrt(2)*omega_c*s+omega_c^2);
Fsd = c2d(Fs,ts,'zoh');
[numfsd,denfsd] = tfdata(Fsd,'v');  % discrete filter
invP = s*(T*s + 1)/K;   % inverse model
Est = Fs * invP         % estimation
Estd = c2d(Est,ts,'zoh');
[numestd,denestd]  =tfdata(Estd,'v'); % discrete model

% Experiment
r      = 60;
r_cyc  = 4;
dist   = 2; % 0 or 2 for dist test
% dist = 0;
Ncyc   = 8;
tfinal = r_cyc*Ncyc;
open_system('pos_pid_mbd_sl_observer_learning')
open_system('pos_pid_mbd_sl_observer_learning/Scope')
sim('pos_pid_mbd_sl_observer_learning')

% EOF of pos_pd_mbd.m

実験結果

まず前回の繰り返し制御を入れていない, PD制御と外乱オブザーバの波形を示します。 最近動かしてなかったとはいえ少し汚い波形でちょっと笑ってしまいました。（ギアもガタガタでグリスを塗りなおした記憶があります）

繰り返し制御を入れた時の波形を示します。 過渡応答の部分を見ると, 繰り返し制御なしの場合は大きなオーバーシュートが発生していることが確認できますが, 繰り返し制御ありの場合は外乱が抑制されオーバーシュートがなくなり, 応答が向上していることがお分かりかと思います。

比較用に, 外乱オブザーバも繰り返し制御もない場合の波形も示します。 この場合は定常的な外乱の影響で指令値に追従できていない状況となります。 再び, 繰り返し制御つきの応答を見てみましょう。 定常外乱が発生してもオーバーシュートなく指令値に追従している波形であることが分かります。

この波形は8試行回した時の波形なので分かりにくいのですが, 応答自体が漸近的に性能向上していく動きになります。(ループをかなり回してみたら結構面白かったです) ループの回数を増やせば増やすほど性能は向上しますが, ゲインや試行回数によっては発振していきます。 時間の都合上, 今回は安定性の考察など設計法を説明することができませんが, 必要ではあると感じています。（次回あればやらせていただきます）

ちなみに僕は学習しすぎると頭が痛くなっちゃうのですが, システムの脳みそであるコントローラもそうなっちゃうのかなと解釈しています。

補足

以下に応答波形のみを載せます。少し時間をおいて取り直してみました。

w/o 外乱オブザーバ, w/o 繰り返し制御 w 外乱オブザーバ, w/o 繰り返し制御 w 外乱オブザーバ, w 繰り返し制御

最終的に1次遅れ系の応答のように落ち着きました。

補足2

動画を撮ってみたのでその様子を載せておきます。

モーション制御 (繰り返し制御)

https://www.youtube.com/watch?v=1SiAKOex_-w&feature=youtu.be

以前外乱オブザーバのみの制御で撮影したものも載せておきます。

モーション制御（外乱オブザーバ）

https://www.youtube.com/watch?v=6yTACA7ZFrU&feature=youtu.be

過渡応答の動きのオーバーシュートが抑えられているかと思います。

感想

”学習”という言葉に関しては少し議論が必要かなと思っていますが, AIのようなものに時間を使いすぎるのもなぁとか, 僕という一個人がそういうワードを使いたくてなんとなく組み込んでみたかっただけなので, あんまり突っかかってこないでいただけると幸いです。(いろいろ言われたら修正する必要はあるとは思いますけど)

そもそもフィードバック制御そのものに上がりすぎたら下げる, 下がりすぎたら上げるといった簡単な学習的な機構があるよね, といった気持ちはあります。

機械学習とか深層学習とか（僕はやったことないんですけど）は, 制御工学の背景知識も含めることによって高度な自動化が達成されるのではないかなと。 大電力を扱っていたり, 人が載ってるわけではないので, 今回の条件なら学習はこんなかんじでいいだろうと思ってます。あいまいさを許容するという意味合いで。

短いですが, 今回はここまでです。今回はあまり準備する時間がなく, 複雑なこともしたくなかったのですが, この制御入れてみたらどうなるだろうかなと思ったらすぐに試せたうえに, 自分の推察を結構簡単に試行錯誤できたりするので, Matlab/Simulink はいいツールじゃないのかなと思っています。

繰り返し制御について, 検索しただけで基礎理論や応用に関する論文がヒットするくらい扱いの多いアルゴリズムのようです。個人的に制御工学 Advent Calender 4日目のほうにも少し書いてみましたので, ご興味がありましたらご覧ください。

繰り返し制御の紹介 https://qiita.com/Manao/items/0bdf01a3fc309a13d721

暗号

暗号文字を解読する

~ のアルファベットを他のアルファベットに暗号化するためのテーブルとして, table[0] ~ table[25] を用意し, を暗号化したときの文字を table[1], ... と格納します。 ~ 以外の文字は table[26] に として格納します。 これを写像としてとらえてみます。 以下にプログラムを示します。

プログラム

#include <stdio.h>

void main(void)
{
    static char table[] = {'Q','W','E','R','T','Y','U','I','O','P',
                           'A','S','D','F','G','H','J','K','L','Z',
                           'X','C','V','B','N','M','?'};
    char *a,*crypto="KSOIDHEPZ";
    int index;

    a=crypto;
    while (*a!='\0'){
        if('A'<=*a && *a<='Z')
            index=*a-'A';
        else
            index=26;
        putchar(table[index]);
        a++;
    }
}

table[] が暗号化テーブルになります。

実行結果

ALGORITHM

組み合わせ（漸化式）

を漸化式を使って求めてみる

ja.wikipedia.org 個の中から 個を選ぶ組み合わせの数を と表現し, 以下のような式で表されます。

は で表される階乗。 はConbination の頭文字から取っています。 この式で計算した場合, 大きな の時に がオーバーフローすることがあります。 int 型なら は なのでintが扱える範囲を超えています。

そこで漸化式を用いて表現することを考えます。漸化式は、各項がそれ以前の項の関数として定まるという意味で数列を再帰的に定める等式です。

例えば は,

という漸化式を利用して表すこともできます。こちらは乗法表示による表現であり、

という純加法表現で表すこともできます。

プログラム

漸化式をプログラムにする場合は、繰り返し又は再帰呼び出しを用いて表現できます。

以下は繰り返しによる乗法表現

#include <iostream>
#include <iomanip>
using namespace std;

long combination(int,int);

int main(){
    int n,r;
    for (n=0;n<=6;n++){
        for (r=0;r<=n;r++)
            cout << setw(3) << n << " C " << r << " = " << setw(2) <<combination(n,r);
        cout << endl;
    }
    return 0;
}

long combination(int n,int r){
    int i;
    long p=1;

    for (i=1;i<=r;i++)
        p=p*(n-i+1)/i;
    return p;
}

次は再帰呼び出しによる加法表現です。

#include <iostream>
#include <iomanip>
using namespace std;

long combination(int,int);

int main(){
    int n,r;
    for (n=0;n<=6;n++){
        for (r=0;r<=n;r++)
            cout << setw(3) << n << " C " << r << " = " << setw(2) <<combination(n,r);
        cout << endl;
    }
return 0;
}

long combination(int n,int r){
    if (r == 0 || r == n)
        return (1);
    else if (r == 1)
        return (n);
    return (combination(n-1,r-1) + combination(n-1,r));
}

実行結果

  0 C 0 =  1
  1 C 0 =  1  1 C 1 =  1
  2 C 0 =  1  2 C 1 =  2  2 C 2 =  1
  3 C 0 =  1  3 C 1 =  3  3 C 2 =  3  3 C 3 =  1
  4 C 0 =  1  4 C 1 =  4  4 C 2 =  6  4 C 3 =  4  4 C 4 =  1
  5 C 0 =  1  5 C 1 =  5  5 C 2 = 10  5 C 3 = 10  5 C 4 =  5  5 C 5 =  1
  6 C 0 =  1  6 C 1 =  6  6 C 2 = 15  6 C 3 = 20  6 C 4 = 15  6 C 5 =  6  6 C 6 =  1

モンテカルロ法

モンテカルロ法を用いてπを求める

ja.wikipedia.org

ある問題を確率（乱数）を用いて解くことをモンテカルロ法といいます。円周率πをこの方法で求めてみます。

~ の一葉実数乱数を2つ発生させ、それらを 、 とします。こうした乱数の組をいくつか発生させると、 の正方形の中に、 で示される点は均一にばらまかれると考えられます。

したがって、正方形の面積と 円の面積の比は、ばらまかれた乱数の数に比例するはずです。

円の中にばらまかれた乱数の数を 、円外にばらまかれた乱数の数を とすると、次の関係が成立します。

プログラム

#include <iostream>
#include <iomanip>
#include <math.h>
#include <stdlib.h>
using namespace std;

#define NUM 330000

double rnd(void);

int main(){
    double x,y,pai;
    int i,in=0;

    for (i=0;i<NUM;i++){
        x=rnd();
        y=rnd();
        if (x*x+y*y<=1)
            in++;
    }
    pai=(double)4*in/NUM;
    cout << "Pi=" << pai << endl;

    return 0;
}
double rnd(void)
{
    return ((double)rand()/32767.1);
}

実行結果

Pi=3.14156