EXCEL2010の散布図テクニック

講義で必要になったので調べたのだが、EXCEL2007から散布図の各要素にラベルがつかなくなった。
これは明らかに不便だ。EXCEL2003まではあったのに。

かつてはこれに以下のようにラベルをつけてどの点が誰なのか分かるようにできた。

講義前日になってこの事実が明らかとなり、
慌ててwebを探し回ったら microsoft のページにわざわざラベル表示専用のマクロが用意されていた。↓
散布図のデータ ポイントにラベルを追加するマクロ
そんなことするくらいなら、機能削除しなければ良かったのに、と文句言いながらシートに実装。
ただ、ほとんどの学生はマクロなど使ったことが無いはずで、
いちいちマクロ機能を使えるようにEXCELを設定しなおして
マクロを自分で呼び出して実行させる、なんてことしたら混乱の極み。
何とかボタン押すだけでラベル付けられるようにならないかと更にwebを渡り歩く。
で、出来上がったのが以下のマクロ。
これは、散布図を作るであろうシートに「ラベル付け」ボタンを付けておき、
このボタンでマクロを起動するようにしておくと、
そのシート内にあるあらゆるグラフを探し出して（散布図に限らず）
ラベルを付けてくれる。
但し、そんなに賢く作ってなくて、散布図データとして囲った
最初の列の一つ前の列をラベルと認識してしまう。だから上の例でいくと、
算数と理科を囲った場合はその前の列のA,B,Cをラベルにつけてくれるが、
理科と社会を囲うと、その前の算数の点をラベルにしてしまう。
とりあえず、サンプルシート(.xlsで保存してあります。)→110529_scatter.xls
分かる人、この先改良してね。

Sub AttachLabelsToPoints()
   'Dimension variables.

    Dim Counter As Integer, ChartName As String, xVals As String

    

    Dim intObj As Integer

    Dim i As Integer

    Dim chart0 As Chart
    'アクティブシートのChartObjects数をカウント

    intObj = ActiveSheet.ChartObjects.Count

   

    'Disable screen updating while the subroutine is run.

    Application.ScreenUpdating = False

    

    'そのシートにあるすべてのグラフをChartとして取得してラベルづけする

    For i = 1 To intObj

    

    Set chart0 = ActiveSheet.ChartObjects(i).Chart

   

    'Store the formula for the first series in "xVals".

    xVals = chart0.SeriesCollection(1).Formula
    'Extract the range for the data from xVals.

    xVals = Mid(xVals, InStr(InStr(xVals, ","), xVals, _

        Mid(Left(xVals, InStr(xVals, "!") - 1), 9)))

    xVals = Left(xVals, InStr(InStr(xVals, "!"), xVals, ",") - 1)

    Do While Left(xVals, 1) = ","

        xVals = Mid(xVals, 2)

    Loop
    'Attach a label to each data point in the chart.

    For Counter = 1 To Range(xVals).Cells.Count

        chart0.SeriesCollection(1).Points(Counter).HasDataLabel = True

        chart0.SeriesCollection(1).Points(Counter).DataLabel.Text = _

            Range(xVals).Cells(Counter, 1).Offset(0, -1).Value

    Next Counter

   

    Next i
End Sub

ガラガラポンのある試み

3年生向けに「統計とコンピュータ」なる講義を毎年やっている。
「確率変数が確率分布に従う」といった概念がなかなか分かってくれなかったり、
そもそも確率変数がピンと来なかったり。
（ガラガラポンとか言っていつも誤魔化しているのだけど）
さてそろそろ、この時期は確率論と統計学を橋渡しする話をする頃。
「アンケート調査」といった統計的な営みをどのように確率へ結びつけていくのか？
標本平均を集めるとどんな現象が起こるのか？
毎年、分かったような分からないような感じのまま推定検定の形式操作に進んでしまい、
いつも振り返って不完全燃焼感が残っていやだった。
そこで今年初めての試みをやってみた。



確率的に30%の人がYes=1といい、残りがNo=0という状況で500×500=25万人の仮想母集団を作り、
各自がその集団の一部（50×50=2500人）にアンケートをとる、
というプログラムを10進BASICで作った。
右図が母集団の一例。Yesと言った人をグレーで塗りつぶしてある。
これが神様の視点だ。

しかし人間にはアンケートによって一部しか測ることができない。
各自はマウスによって任意の場所に移動して仮想アンケートを取る。（左クリック！）
するとこの人間の目に見える2500人中に何人Yesといったかをプログラムが数えるようにしてある。
こういった実験を各自50回ずつ行ってその結果をあるwebページへ送信してもらい、
web上で集計しヒストグラムにする。
利用したサイトは「Mr.アンケート」。無料で手軽に使える。



結果は右図。3クラスやったが、
A組B組はきれいになった（素直だ）。
そして結果を見れば分かるように、C組には遊ばれた（泣）
始めはギザギザだったヒストグラムも、
時間とともに調査データが集まりだすと
思ったよりもきれいに正規分布の形状をなしてきたのは嬉しい。
ま、二項分布B(2500,0.3)に従う確率変数データを集めたのだから、きれいになりやすいのも尤もなんだが。


仮想ではあったが各自が「アンケート調査」し、
それを一箇所に集約した結果、正規分布が現れてくる様子を目の当たりにすることで、
幾らかは統計と確率論が学生の中で結びついてくれるといい。
そう、「アンケートに一人が答える」＝「ガラガラポンを一回振る」
となってくれているだろうか？


アンケートシミュレーションで使ったBASICプログラム⇒[100609enquete.bas]

【使い方】
実行すると仮想25万人のYes/Noの様子が打ち出され、
その後グレーの「無知のヴェール」で場が覆われます。
一部50×50=2500人分だけ窓が開いていてそこに入っている人たちに
左クリックでアンケートをとり、Yesの人数が表示されます。
場所を移動しつつ左クリックすることで次々に2500人アンケート結果が表示されます。
飽きたら右クリックで終了。

REM 【統計とコンピュータ】

REM 内閣支持率の電話アンケートシミュレーション

REM
LET rm=.3            !理論的な支持率

LET sc=5             !一マス=一人の大きさ

LET L=500            !window一辺のマスの数(人数)

LET sz=L*sc          !windowの一辺の長さ

LET dd=sc-1

LET pop=L*L          !全人口

LET dnum=pop*rm      !理論的にいるはずの全支持者数

LET asz=.1           !アンケートサイズ

LET apop=(L*asz)^2   !アンケート人数
OPTION BASE 0

DIM dots(L,L)
PRINT "全人口は";L;"×";L;"=";pop;"人"

PRINT "理論的な支持率=";rm

PRINT "1回の調査人数=";apop

PRINT "理論的には調査で平均";apop*rm;"人支持者がいるはず"

PRINT

PRINT "以下，実際のアンケート調査結果--------------"

PRINT
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

! 支持者のランダム配置

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SET WINDOW 0,sz,0,sz
SET LINE COLOR 4

PLOT LINES: 0,0;sz,0;sz,sz;0,sz;0,0

SET LINE COLOR 8

SET AREA COLOR 1

RANDOMIZE

FOR cx=0 TO L-1

   FOR cy=0 TO L-1

      IF RND<rm THEN

         DRAW box WITH SHIFT(cx*sc,cy*sc)

         LET dots(cx,cy)=1

      END IF

   NEXT cy

NEXT cx
gsave "dots.bmp"
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

! 電話アンケート（調査人数=(L/10)^2）

! マウスで適当に場所を移動してクリックすると

! 囲まれた範囲の支持者の数が表示される

! 右クリックで調査終了

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SET LINE COLOR 4

SET AREA COLOR 8

SET COLOR MIX(8) .9,.9,.9

LET cc=1

DO UNTIL(right=1)

   mouse poll mx, my, left, right

   IF (mx<>mx0 OR my<>my0) THEN

      IF mx>.9 THEN LET mx=.9

      IF my>.9 THEN LET my=.9

      LET ax=INT(mx*L)/L

      LET ay=INT(my*L)/L

      SET DRAW mode hidden

      gload "dots.bmp"

      PLOT LINES: ax,ay; ax+asz,ay; ax+asz,ay+asz; ax,ay+asz; ax,ay

      paint ax/2,ay/2

      SET DRAW mode explicit

      LET mx0=mx

      LET my0=my

   END IF

   IF left=1 THEN

      LET sx=ax*L

      LET sy=ay*L

      LET count =0

      FOR i=0 TO L*asz-1

         FOR j=0 TO L*asz-1

            IF dots(i+sx,j+sy)=1 THEN LET count=count+1

         NEXT j

      NEXT i

      PRINT USING "###回目 ####":cc,count

      LET cc=cc+1

      WAIT DELAY .5

   END IF

LOOP
gload "dots.bmp"
PICTURE box

   PLOT LINES: 0,0; 0,dd; dd,dd; dd,0; 0,0

END PICTURE

END

素朴なＧＡ

今年のゼミ生の卒論の一つで遺伝的アルゴリズム（ＧＡ）についての
素朴な考察を行った。
その際、実際にどれほど有効なものなのか見ようと、
BASICでアルゴリズムを組んでみた。
遺伝子長は4で長さ2に固定した2点交叉
（つまり連続した2遺伝子を両親の間で交換する）
のみで、突然変異などはない。
各染色体が親として選ばれる確率は適合度関数に比例させてある。
下の図は区間[0,1]を16等分して各区間を長さ4の染色体で対応させ
[0,1]区間上の正値関数をそのまま適合度関数にとった。
例えば

f(x)=30x(x-0.3)(x-0.6)(x-0.8)(1-x)+0.5

は染色体が対応している座標7/16=[0111]と15/16=[1111]で
きわめて近い値で極大値となっているが
僅かにf(7/16)>f(15/16)となっている。
(f(7/16)/f(15/16)=1.01412...)
実行結果の以下の図を見ると初めは染色体[1111]も
それなりの確率を持っているが、
400代ほど世代交代をさせると[0111]に確率が集中してくる。


う〜ん、面白い。
どうやら、引き伸ばしと折りたたみの仕組みが
働いているようだ。

10進BASICプログラム→[ga00.bas]

REM [GA00.bas]

REM Genetic Argorithm

REM 遺伝子長=4, 交叉長=2
!DEF f(x)=4*x*(1-x)      ! 適合度関数

!DEF f(x)=30*x*(x-.5)*(x-.5202314098)*(1-x)+.2      ! 適合度関数

DEF f(x)=30*x*(x-.3)*(x-.6)*(x-.8)*(1-x)+.5      ! 適合度関数

!DEF f(x)=EXP(-(x-.5)^2/.002)      ! 適合度関数

LET  win=.02            ! windowの余白

LET L=4                 ! 遺伝子長

LET pt=2^L              ! 遺伝子パターン数

LET dx=1/pt

OPTION BASE 0

DIM g$(pt)              ! 遺伝子名(Binary)

DIM p(pt),q(pt),a(pt),aaoo(pt),ooaa(pt),aooa(pt),oaao(pt)
REM 【初期確率、適合度値】

FOR k=0 TO 15

   LET q(k)=1/pt

   LET p(k)=0

   LET a(k)=f(k*dx)

NEXT k
REM 【交叉パターンの設定】
DIM g1(4,4),g2(4,4),g3(4,4),g4(4,4)

FOR j=0 TO 3

   FOR i=0 TO 3

      LET g1(j,i)=j*4+i

      LET g2(j,i)=j+i*4

      LET g3(j,i)=INT(j/2)*8+MOD(j,2)+i*2

      LET g4(j,i)=j*2+INT(i/2)*8+MOD(i,2)

   NEXT i

NEXT j

FOR k=0 TO 15

   LET aaoo(k)=INT(k/4)

   LET ooaa(k)=MOD(k,4)

   LET aooa(k)=INT(k/8)*2+MOD(k,2)

   LET oaao(k)=INT(MOD(k,8)/2)

NEXT k
REM 【適合度関数の表示】
SET VIEWPORT win, 1-win, .5+win, 1-win

SET WINDOW -.05,1.05,-.1,1.1

SET TEXT COLOR 8

PLOT TEXT ,AT 0,1:"適合度関数"

SET TEXT font "Times New Roman Italic",7

DRAW grid0(dx,.5)
SET LINE COLOR 2

SET LINE width 2

FOR x=0 TO 1 STEP dx

   PLOT LINES : x,f(x);

NEXT x

FOR k=0 TO pt-1

   FOR i=L TO 1 STEP -1

      LET g$(k)=g$(k)&STR$(INT(MOD(k,2^i)/2^(i-1)))

   NEXT i

   PLOT TEXT ,AT k*dx,-.05*(MOD(k,2)+1):g$(k)

NEXT k
REM 【遺伝子存在率のヒストグラム】
SET VIEWPORT win, 1-win, win, .5-win

SET WINDOW -.05,1.05,-.1,1.1

SET LINE width 1

SET TEXT font "ＭＳ ゴシック",11

PLOT TEXT ,AT 0,1:"各遺伝子パターンの存在確率"

SET TEXT font "Times New Roman Italic",7

DRAW grid0(dx,.5)

PLOT TEXT ,AT -.03,1:"1.0"

PLOT TEXT ,AT -.03,.5:"0.5"

FOR k=0 TO pt-1

   PLOT TEXT ,AT k*dx,-.05*(MOD(k,2)+1):g$(k)

NEXT k
INPUT PROMPT "世代数=":gen

LET gstep=INT(gen/10)
REM 【交叉計算とヒストグラム表示】
SET LINE COLOR 10
FOR t=1 TO gen

   LET S=0

   FOR k=0 TO 15

      LET S=S+a(k)*q(k)

   NEXT k

   FOR k=0 TO 15

      LET p(k)=a(k)*q(k)/S

   NEXT k

   LET p1=0

   LET p2=0

   LET p3=0

   LET p4=0

   FOR k=0 TO 15

      LET p1=p(g1(aaoo(k),0))+p(g1(aaoo(k),1))+p(g1(aaoo(k),2))+p(g1(aaoo(k),3))

      LET p2=p(g2(ooaa(k),0))+p(g2(ooaa(k),1))+p(g2(ooaa(k),2))+p(g2(ooaa(k),3))

      LET p3=p(g3(aooa(k),0))+p(g3(aooa(k),1))+p(g3(aooa(k),2))+p(g3(aooa(k),3))

      LET p4=p(g4(oaao(k),0))+p(g4(oaao(k),1))+p(g4(oaao(k),2))+p(g4(oaao(k),3))

      LET q(k)=(p1*p2+p3*p4)/2

   NEXT k

    

   IF MOD(t,gstep)=0 OR t=gen THEN

      SET COLOR MIX(10) (1-t/gen)/2,(1-t/gen)/2+.5,t/gen    

      FOR k=0 TO 15 

         PLOT LINES : k*dx,q(k);

      NEXT k

      PLOT lines

   END IF

    

NEXT t

END

プリアンブルの不思議

LaTeXでjpg画像などを取り込んでプリントを作る機会が多いのだが、
webにもアップする為pdfにまで変換する。
が、時々dvioutまでは画像として出るのに、
pdfにはまり込んでないことがある。
前回とほとんど同じファイルから出発したのに、
一方はpdfまで画像が出せ、一方は落ちてしまう。


随分長い間原因が分からなかったのだが、
今日実験していて、ちょっとだけ分かったので記録しておく。
どうやらプリアンブルにパッケージを読み込む順番に
問題があったようだ。


失敗していたのは以下の状態。

\usepackage{tabularx,ascmac}
\usepackage[dviout,dvipdfm]{graphicx}
\usepackage{amsmath}
\usepackage{float}

うまくいっていたのが以下の状態。

\usepackage[dvipdfm]{graphicx}
\usepackage{amsmath}
\usepackage{ascmac}

どれかのパッケージ同士がかち合ってしまっているんだろうか。
ま、とにかくうまくいかなかったら、
読み込み順に気をつけてみよう。

フーリエでモアレ

フーリエ級数と周期性ノイズについてまとめている学生のネタ。
2次元のうなりとしてのモアレを
わざわざフーリエ変換したデータを傷つけることによって
モアレを擬似的に起こすシミュレーションをした。
以下がそのプログラム。
やっつけで作っているので、相当汚い。


上段が元画像とそのフーリエ変換、
下段がフーリエ変換データに傷を付けたものと
その逆フーリエ変換。

モアレノイズがどう消せるのかの
一つの考察。


こんな感じ↓





10進BASICプログラム→[moare.bas]

REM [FFT2D_moare.bas]
REM 2D FFT/IFFT でモアレを合成する
REM FFT/IFFTは外部関数として使う

OPTION ARITHMETIC complex
SET TEXT FONT "",6
SET AXIS COLOR 0

DECLARE EXTERNAL SUB FFT    ! FFT
DECLARE EXTERNAL SUB IFFT   ! IFFT

! 配列サイズ dn は必ず2^nで
LET dn = 256            ! データ総数
LET  dn2 = dn/2
LET  win=.02            ! windowの余白

LET  dt = 1   ! サンプリング間隔
LET  ct= 10^0  !元画像のコントラスト
LET  mct= 10^(.5)  !モアレ画像のコントラスト
LET  fct= 10^(-.5)  !フーリエスペクトルのコントラスト

!極座標で(r,θ)=±(mr,mth)の周波数のところにmstの強さで傷を付ける
LET  mst=100  !モアレ傷の強さ
LET  mth=30   !モアレの角度
LET  mr=.5 !モアレの細かさ

DEF f(x,y)=(1-SGN(x^2+y^2-.5)-.1)+10*(1-SGN(2*(x-.2)^2+(y-.2)^2-.01)-.1)
　+10*(1-SGN(2*(x+.2)^2+(y-.2)^2-.01)-.1)
　+10*(MAX(-SGN(x^2-y-.3)+SGN(x^2/2-y-.2)-.1,0))

SET POINT COLOR 10
SET POINT STYLE 1
SET bitmap SIZE dn*2*(1+win*2),dn*2*(1+win*2)

REM 【元の画像の表示】

SET VIEWPORT win,.5-win , .5+win,1-win
DIM bmp(dn,dn)          ! 画像用配列 dn×dn

SET WINDOW -dn2, dn2, -dn2, dn2
FOR y=-dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET  g=f(x/dn2,y/dn2)
      LET bmp(x+dn2+1,y+dn2+1) = g
      LET  g=ATN(g*ct)/PI+.5
      SET COLOR MIX(10) g,g,g 
      PLOT POINTS: x,y
   NEXT x
NEXT y

CALL FFT(bmp)

REM 【周波数スペクトラム】

SET VIEWPORT .5+win,1-win , .5+win,1-win
SET WINDOW -dn2, dn2, -dn2, dn2
FOR y= -dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET g=ATN(ABS(bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1))*fct)/PI*2
      SET COLOR mix(10) g,g,g
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

REM 【モアレ傷を付ける】

SET VIEWPORT .5+win,1-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
LET  mx=mr*COS(mth/180*PI)
LET  my=mr*SIN(mth/180*PI)

FOR y= -dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      IF (ABS(x-dn2*mx)<1 AND ABS(y-dn2*my)<1) THEN
         LET  bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1)=mst
      END IF
      IF (ABS(x+dn2*mx)<1 AND ABS(y+dn2*my)<1) THEN
         LET  bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1)=mst
      END IF
      LET g=ATN(ABS(bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1))*fct)/PI*2
      SET COLOR mix(10) g,g,g
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

CALL IFFT(bmp)

REM 【モアレ画像の表示】
SET VIEWPORT win,.5-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2

FOR y= 1 TO dn
   FOR x = 1 TO dn
      LET g=ATN(re(bmp(x,y)+.1)*mct)/PI/2+.5
      SET COLOR mix(10) g,g,g
      PLOT POINTS: x-dn2,y-dn2
   NEXT x
NEXT y

END

REM  [2D FFT/IFFT]
REM   元データ x(,) → FFT/IFFT → 変換後データ x(,)
REM   ユニタリ変換になるようFFT/IFFT両方を1/N
REM   Fk = 1/N * Σ[n=1...N] (fn * W^( (k-1)*(n-1)) )   :FFT
REM   fk = 1/N * Σ[n=1...N] (Fn * W^(-(k-1)*(n-1)) )  :IFFT
REM   W = exp(-��(-1) * 2π/N)
REM   配列添字は 1 から出発

REM 【2D FFT】
REM θ = -2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB FFT(x(,))
!OPTION ARITHMETIC NATIVE
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = -2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM 【2D IFFT】
REM θ = 2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB IFFT(x(,))
!OPTION ARITHMETIC NATIVE
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = 2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM [バタフライ演算]
REM 配列サイズを半分にしながら再帰呼び出しする
REM wk= exp( (k-1)θ )
REM x(k)    --- *1  -- x0(k)
REM         × 
REM x(k+N/2)--- *wk -- x1(k)

EXTERNAL SUB FFTMAIN(x(), th)
!OPTION ARITHMETIC NATIVE
OPTION ARITHMETIC complex
LET sx = SIZE(x)
IF sx>1 THEN
   LET sx2 = sx/2
   DIM x0(sx2), x1(sx2)
   FOR k=1 TO sx2
      LET x0(k) = x(k) + x(k+sx2)
      LET wk = EXP( complex(0, th * (k-1)) )
      LET x1(k) = wk * ( x(k) - x(k+sx2) )
   NEXT k
   LET th0 = 2 * th
   CALL FFTMAIN(x0, th0)
   LET th1 = 2 * th
   CALL FFTMAIN(x1, th1)
   FOR k=1 TO sx2
      LET x(2*k-1) = x0(k)
      LET x(2*k) = x1(k)
   NEXT k
END IF
END SUB