ニュートンの第2法則


第一法則では、質点の速度ベクトルに何らかの変化が認めたならば、力が作用したと解釈する。これは力の定性的な定義として考えることができる。
第二法則では、運動量の変化を力としてより定量的に定義を行っている。




この流れからわかるよに、第2法則は力の定義という意味合いだけでなく、既知の力のもとで質点の運動を予測する方程式としての意味も持ってくるのである。

ちなみにm=constの前提は、ロッケトが燃料を消費して飛び立つ時などは成り立たちません。l

%text.tex

\documentclass{jsarticle}
\begin{document}
Momentum
\begin{eqnarray}
\mbox{\boldmath$p$} \equiv m \mbox{\boldmath$v$}
\end{eqnarray}

Equation of motion
\begin{eqnarray}
\frac{d \mbox{\boldmath$p$}}{dx}=\mbox{\boldmath$F$}\\
m \frac{d \mbox{\boldmath$v$}}{dt}+ \frac{dm}{dt}\mbox{\boldmath$v$}=\mbox{\boldmath$F$}
\end{eqnarray}
We assume that Mass dose not depend on Time.  (m = const)
\begin{eqnarray}
m \frac{d \mbox{\boldmath$v$}}{dt}=\mbox{\boldmath$F$}\\
m \frac{d^2 \mbox{\boldmath$r$}(t)}{dt^2}=\mbox{\boldmath$F$}
\end{eqnarray}
The initial conditions are needed to solve this equation.
\begin{eqnarray}
\mbox{\boldmath$r$} = \mbox{\boldmath$r_{0}$},\frac{d\mbox{\boldmath$r$}}{dt} = \mbox{\boldmath$v_{0}$}(t = t_{0})
\end{eqnarray}
or
\begin{eqnarray}
\mbox{\boldmath$r$} = \mbox{\boldmath$r_{1}$}(t = t_{1}), \mbox{\boldmath$r$} = \mbox{\boldmath$r_{2}$}(t = t_{2}),
\end{eqnarray}
Given the place and velocity, we solve the equation uniquely.
These are the law of cause and effect in classical mechanics.


\pagestyle{empty}
\end{document}

ニュートンの第2法則

ニュートンの3法則では、第一法則において「力」を質点の慣性に基づく運動に影響を与えるものとして定性的に定義されています。

そして、第2法則では力を質点の運動量の変化として定量的に定義されています。


ニュートンの運動の3法則

藤原邦男先生の「物理学序論としての力学」に記載されている運動の3法則をメモしておく。

法則1
すべての質点はそれに加えられた力によってその状態が変化させられない限り、静止あるいは1直線上の等速運動の状態を続ける。

法則2
質点の運動量(=質量×速度)の変化は、加えられた力の方向にそって起こり、かつ、微小時間内における運動量の単位時間あたりの変化の大きさは、加えられた力の大きさに等しい。

法則3
すべての作用に大して、等しく、かつ繁体向きの反作用が常に存在する。すなわち、お互いに働き合う2つの質点の相互作用は常に相等しく、かつ反対方向へと向かう。

英語だと以下。

$ vim text.tex


\documentclass{jsarticle}
\begin{document}
Newton's Three Laws of Motion
\begin{enumerate}
\item Every object in a state of uniform motion tends to remain in that state of motion unless an external force is applied to it.
Probability density function
\item The relationship between an object's mass m, its acceleration a, and the applied force F is F = ma. Acceleration and force are vectors (as indicated by their symbols being displayed in slant bold font); in this law the direction of the force vector is the same as the direction of the acceleration vector.
\item For every action there is an equal and opposite reaction.
\end{enumerate}

\pagestyle{empty}
\end{document}

$ platex text.tex
$ dvipng -D 500 -T tight text.tex



統計学の基本的な関数の定義

$ text.tex

% text.tex
\documentclass{jarticle}

\begin{document}

Probability density function
\begin{eqnarray}
P(a \leq X \leq b) = \int_a^b f(x)dx \\
f(x) \geq 0 \ \ (- \infty \leq x \leq \infty) \\
\int_{-\infty}^{\infty} f(x)dx = 1
\end{eqnarray}

Cumulative distribution function
\begin{eqnarray}
F(x) = \int_{-\infty}^x(fu)du \\
\frac{d}{dx}F(x) = f(x)
\end{eqnarray}

Expected value
\begin{eqnarray}
E(x) = \int_{-\infty}^{\infty}xf(x)dx = \mu \\ 
\end{eqnarray}

Variance
\begin{eqnarray}
V(x) = E \{ (X - \mu ) \}^2 = \int_{-\infty}^{\infty}(x - \mu)^2f(x)dx = \sigma^2 
\end{eqnarray}

\pagestyle{empty}
\end{document}



$ platex text.tex
$ dvipng -D 1000 -T tight text.dvi







指数分布

故障率が一定のシステムで次の故障が起きるまでの時間Xや次の災害が起きるまでの時間Xなどのまったく独立に生じる現象は、その待ち時間が指数分布に従うことが知られています(2)。


待ち時間の累積分布関数は(2)のようになる。

これらのことは、指数分布に従う現象は、近い将来に起こっても不思議ではないことを意味する。
航空機の墜落のような希少事象が連続することがこれに該当する。

$ vim exponential_distribution.tex



%exponential_distribution.tex


\documentclass{jsarticle}
\begin{document}

Probability density function of exponential distribution
\begin{equation}
f(t)= \left \{
\begin{array}{l}
f(t) = \lambda \mathrm{e}^{-\lambda t}\ (t \geq 0)\\
0 \ ( t < 0 )
\end{array}
\right.
\end{equation}

Cumulative distribution function of exponential distributiont   $ (t\geq 0) $
\begin{equation}
F(t) = P(X \leq t) = \int_{-\infty}^t f(u) du = 1 - \lambda \mathrm{e}^{-\lambda t}
\end{equation}

\end{document}

$ platex exponential_distribution.tex

$ dvipng -T tight exponential_distribution.dvi






f(t) , F(t)をRを用いて描写してみます。
$ vim command.in

png("120816_figure.png")
par(mfrow=c(1,2))
curve(exp(-1 * rambda * x),0,10,main="f(t), rambda = 1")
curve(1 - exp( -1 * rambda * x) ,0,10,main="F(t), rambda = 1")
dev.off()
history()

$ R
source("command.in")


scpコマンドでファイルを転送

ファイルの内容を暗号化して、ftpよりも安全に他のコンピュータに送るときは、scpを使用します。


$ scp hoge.file username@ipaddress:path

運動方程式を題材にしてplatexによる数式のフィギアを作る。


○基本形は下

\documentclass{jsarticle}
\begin{document}
%コメントをここに入れる
\pagestyle{empty}
\end{document}


○改行して数式を入れる


\documentclass{jsarticle}
\begin{document}
\[%ここに数式 \]
\pagestyle{empty}
\end{document}


○太字斜体でベクトルを表記する
\mbox{\boldmath $a$}


○時間微分
\dot a




○二階の時間微分
\ddot a



○コンパイル
$ platex hoge.tex


○dviをpngに変換
$ dvipng -T tight hoge.dvi

それでは、いよいよニュートンの運動方程式をtexで作成してみます。

$ vim newton.tex


% newton.tex
\documentclass{jsarticle}
\begin{document}
\[ m\frac{d^2\mbox{\boldmath $r$}}{dt^2}= \mbox{\boldmath $F$} \]
\pagestyle{empty}
\end{document}



$ platex newton.tex

$ dvipng -T tight newton.dvi




















platexのインストールおよび簡単な使い方

$ vim platexinstall.sh


# platexinstall.sh
sudo apt-get install -y texlive texlive-math-extra ptex-bin xdvik-ja

sudo apt-get install -y dvipdfmx cmap-adobe-japan1 ptex-jisfonts okumura-clsfiles jmpost jbibtex-base jbibtex-bin mendexk

sudo jisftconfig add


$ sudo chmod +x latexinstall.sh

$ sudo ./latexinstall.sh


#テストスクリプトをかく
$ vim hello.tex


\documentclass{jsarticle}
\begin{document}
Hello, world!
\pagestyle{empty}
\end{document}


#コンパイル。dviファイルが生成する。
$ platex hello.tex

#dvipngをインストール
$ sudo apt-get install dvipng

#dviファイルをpngに変換(余白はタイトに切り詰める)
$ dvipng -T tight helllo.dvi




\pagestyle{empty}を入れておかないと、デフォルトでページ番号が入ってしまい、pngにしたときに、文頭からページ番号までを含んだ縦長のファイルが生成してしまうので注意が必要である。







標準偏差一つ分、二つ分。。。

小児の発達などで登場する-2SDという概念を考察します。

SD = Standard Deviation(標準偏差)

ここで、標準偏差1、平均値0の標準正規分布の確率密度曲線を考えます。

# -SD - SDの区間の確率を計算
> integrate(dnorm,-1,1)
0.6826895 with absolute error < 7.6e-15
# -2SD - 2SDの区間の確率を計算
> integrate(dnorm,-2,2)
0.9544997 with absolute error < 1.8e-11

以上の結果から、母集団の-2SDから2SD以内というのは、母集団のおおよそ95%以内に入るということだと理解できます。

参考に図を描写しておきます。

png("120812_curve.png")
curve(dnorm(x),xlim=c(-5,5))
abline(v=-2,col="green")
abline(v=-2,col="green")
abline(v=-2,col="magenta")
abline(v=2,col="magenta")
dev.off()



ヒストグラムの階級を細かくして正規分布の密度関数を近似的に描写する


#標準正規分布から乱数を取得
dat <- rnorm(1000000)
#レンジを調べる
range(dat)
png("120812_hist.png")
par(mfrow=c(2,2))
hist(dat, breaks=seq(-10,10,10))
hist(dat, breaks=seq(-10,10,1))
hist(dat, breaks=seq(-10,10,0.1))
hist(dat, breaks=seq(-10,10,0.01))
dev.off()

dat <- rnorm(1000000)
png("120812_hist_2.png")
par(mfrow=c(1,2))
hist(dat,breaks=seq(-6,6,0.1))
curve(length(dat)*dnorm(x))
dev.off()









カイ二乗検定で独立性の検定を行う

いま、20人の大学生がいます。
それぞれに数学と物理学が好きか嫌いかを二択で答えてもらうとします。
結果が以下のようになったときに、数学と物理の好き嫌いの傾向は独立かどうかを検定します。
$ R
> 観測度数 <- c(10, 2 , 4, 4) 


> 期待度数 <- c(8.4, 3.6, 5.6, 2.4)
> カイ二乗要素 <- (観測度数-期待度数)^2/期待度数
> カイ二乗要素
[1] 0.3047619 0.7111111 0.4571429 1.0666667
> カイ二乗値 <- sum(カイ二乗要素)
> カイ二乗値
[1] 2.539683


> pchisq(カイ二乗値, 1, lower.tail=FALSE)
[1] 0.1110171

p > 0.05であるため、このアンケートの結果は偶然でも起こりえる程度のものであるという結論となります。

ちなみに、df= 1のカイ二乗分布は以下のようになります。

png("120807_chi_squere.png")
curve(dchisq(x,1), 0, 30)
abline(v=qchisq(0.05,1,lower.tail=FALSE), col = "red")
abline(v=カイ二乗値)
dev.off()




Mac OS XにEMBOSSをインストール

UNIX環境における配列解析の定番ソフト「EMBOSS」をビルドしてみます。

#ソースコードのダウンロード
ftp ftp://emboss.open-bio.org
> cd pub
> cd EMBOSS
> get EMBOSS-6.5.7.tar.gz
> quit

#解凍
$ tar xzvf EMBOSS-6.5.7.tar.gz
$ cd EMBOSS-6.5.7

#configureスクリプトを実行
$ ./configure
#コンパイル
$ make
#インストール
$ sudo make install

#試しに、embleからxlrhodopの塩基配列(Xenopus laevis rhodopsin mRNA, complete cds)を取得する。

$ seqret
Read and write (return) sequences
Input (gapped) sequence(s): embl:xlrhodop
output sequence(s) [l07770.fasta]:
kappa-no-MacBook:120806 kappa$ ls
l07770.fasta

>L07770 L07770.1 Xenopus laevis rhodopsin mRNA, complete cds.
ggtagaacagcttcagttgggatcacaggcttctagggatcctttgggcaaaaaagaaac
acagaaggcattctttctatacaagaaaggactttatagagctgctaccatgaacggaac
agaaggtccaaatttttatgtccccatgtccaacaaaactggggtggtacgaagcccatt
cgattaccctcagtattacttagcagagccatggcaatattcagcactggctgcttacat
gttcctgctcatcctgcttgggttaccaatcaacttcatgaccttgtttgttaccatcca
gcacaagaaactcagaacacccctaaactacatcctgctgaacctggtatttgccaatca
cttcatggtcctgtgtgggttcacggtgacaatgtacacctcaatgcacggctacttcat
ctttggccaaactggttgctacattgaaggcttctttgctacacttggtggtgaagtggc
cctctggtcactggtagtattggccgttgaaagatatatggtggtctgcaagcccatggc
caacttccgattcggggagaaccatgctattatgggtgtagccttcacatggatcatggc
tttgtcttgtgctgctcctcctctcttcggatggtccagatacatcccagagggaatgca
atgctcatgcggagtagactactacacactgaagcctgaggtcaacaatgaatcctttgt
tatctacatgttcattgtccacttcaccattcccctgattgtcatcttcttctgctatgg
tcgcctgctctgcactgtcaaagaggctgcagcccagcaacaggaatctgctaccaccca
gaaggctgagaaagaggtcaccagaatggttgttatcatggtcgttttcttcctgatctg
ttgggtgccctatgcctatgtggcattctacatcttcacccaccagggctctaactttgg
cccagtcttcatgaccgtcccagctttctttgccaagagctctgctatctacaatcctgt
catctacattgtcttgaacaaacagttccgtaactgcttgatcaccaccctgtgctgtgg
aaagaatccattcggtgatgaagatggctcctctgcagccacctccaagacagaagcttc
ttctgtctcttccagccaggtgtctcctgcataagagcttcaccagggctgtctcagggt
ccgctgcctcacacaattcccatcacttaagccctgtctacttgttgcgaaggcaaagaa
ttccacagttttaatatttacccccattctgcccaaccttggacactgtaagagctgacc
ccattactgctgggaaggcccaagctttgttgcattctgatgtgatcctttcagcagaaa
atgggtggattcaatgaatttcaccaaggctgtacataacaataacattagtctgaaggc
acctcccacccagagaatgcaacacttatttatctctgtcttttcttgacatattgatgc
tgcttctattcatggtcactaacaaaaagtcccattttacaatgcaactgaaagtaatgt
atttttgtaatataataacatatttcatgcaatctcctctgcttattggcaaggtctgat
atagtgaggatagacagccagaccccttgcattaaaatcctgtattaaaaatttctttgc
aagt


サイコロの出る目の度数分布


各数字が1/6の確率で得られる理想的なサイコロがあるとします。
ここで、このサイコロを繰り返し振って、出た目の度数分布を作成します。


必要な関数(ceiling, runif)の仕様を勉強。

$ R
> ? ceiling
Description:
     ‘ceiling’ takes a single numeric argument ‘x’ and returns a
     numeric vector containing the smallest integers not less than the
     corresponding elements of ‘x’.
与えたnumeric vector(要素は小数でも可)の各要素の小数点を切り捨てた整数値を返します。

>?runif
Description:
     These functions provide information about the uniform distribution
     on the interval from ‘min’ to ‘max’.
与えられた最小値から最大値の間で、一様分布に従う値を返すとのこと。

$ vim command.in

png("120805_dice.png")
par(mfrow=c(2,2))
dice1 <- ceiling(runif(n=6, min = 0, max = 6)) 
dice2 <- ceiling(runif(n=60, min = 0, max = 6)) 
dice3 <- ceiling(runif(n=600, min = 0, max = 6)) 
dice4 <- ceiling(runif(n=6000, min = 0, max = 6)) 
hist(dice1,main="n=6",breaks = seq(0,6,1))    
hist(dice2,main="n=60",breaks = seq(0,6,1))    
hist(dice3,main="n=600",breaks = seq(0,6,1))    
hist(dice4,main="n=6000",breaks = seq(0,6,1))   
dev.off()

$ R
> source("commnad.in")

n数を多くすれば母集団の分布を反映したものが得られます。

検出力と例数設計

Rには、t検定の検出力を計算するための関数が用意されています。

それが、power.t.test()です。

今回はこの関数を用いて、検出力と例数設計について検証してみたいと思います。

先日、身長平均身長175cm(±5cm)と他の集団との身長の平均値に関しての比較を扱いました。その時、有意差を検出するためのn数について直感的な検討を行いました。

今回は、パワーテストなので、より定量的に議論したいと思います。

以下のスクリプトは、1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差 2)想定される平均値の差、標準偏差、α危険率、検出力をある値に想定したときに要請されるn数 3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値 について着目しています。

今回はベースに「平均身長175cm(±5cm)の集団と平均身長180cm(±5cm)の集団の平均値の比較を意識していて、サンプル数は10、p値は0.05を想定している」状況があります。


$ vim command.in


#以下の条件を満たすのに必要なn数を求める。
n_null = power.t.test(delta=5,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件で有意差が検出可能な母平均の差を求める。
delta_null = power.t.test(n=10,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件での検出力を求める。
power_null = power.t.test(delta=5,n=10,sd=5,sig.level=0.05,type=c("two.sample"),alternative="two.sided")


$ R
> source("command.in")

#特定の項目を取り出す方法を調べる。

> names(delta_null)
[1] "n"           "delta"       "sd"          "sig.level"   "power"    
[6] "alternative" "note"        "method"    


#以下、1)-3)のケースの計算結果を表示


> delta_null$delta
[1] 6.624728
> n_null$n
[1] 16.71477
> power_null$power
[1] 0.5619846



1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差


> delta_null$delta
[1] 6.624728
サンプル数は10、p値は0.05、検出力0.8、標準偏差5を想定しているで検出可能な平均値の差は6.6cmであるということです。


2)標準偏差、α危険率、検出力をある値に想定したときに要請されるn数


> n_null$n
[1] 16.71477
平均値の差が5、p値は0.05、検出力0.8、標準偏差5の条件を満足させるために必要なn数は16.7である。
#実験計画を立てる上で特に重要である。



3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値


> power_null$power
[1] 0.5619846
n数が10、平均値の差が5、p値は0.05、標準偏差5の条件の時の検出力は0.56である。


2種類の過誤を明確に意識し、さらには母集団の平均値の差を実験によりある程度想定しながら、サンプル数を決定したいと思いました。

・第一種過誤(α過誤、偽陽性)は、帰無仮説が実際には真であるのに棄却してしまう過誤
・第二種過誤(β過誤、偽陰性)は、対立仮説が実際には真であるのに帰無仮説を採用してしまう過誤


最後に実際に実験を行うことを意識して
「想定する母集団(標準偏差5)の差を1から0.5cmずつ10cmまで増加させていったときに、p値は0.05、検出力0.8の条件を満足させるために必要なn数を2次元プロットに描写してみます。

$ vim command_2.in

n <- c(1:100)
for(i in 1:100)
{
n[i] <- (power.t.test(delta=0.1*i,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided"))$n
}

delta <- c(1:100)*0.1
png("120801_delta_n_plot.png")
plot(delta,n)
dev.off()

$ fg
> source("command_2.in")
このシミュレーション結果は非常に深い意味があります。
平均値の差が小さい時(2ぐらいまで)は、非常に多くのサンプル数が必要となるのに対して、それよりも大きい平均値の差であるときには、サンプル数がほぼ一定(5程度)で済んでしまうというわけです。

deltaの値を0.11から10までの間に関してプロットを書かせるとよりわかりやすいです。

> plot(delta[11:100],n[11:100])



実験を計画するときはこのような事実を念頭に置いてサンプル数を設定したいものです。