読者です 読者をやめる 読者になる 読者になる

ピエール瀧になりたい

備忘録として

統計の勉強-共分散

統計の問題集をちまちま解いているのですが共分散の所でちょっと詰まった。

問題としては独立な二つの正規分布(norm1, norm2)があり、それらの分散は何れも1。

このとき、norm1と(norm1+norm2)/√2の共分散を求めよ、という感じ。

回答みても直感的に理解できなかったのでrで書いてみる。

 

#norm1 <- rnorm(10000, mean = 100, sd = 1) 
#norm2 <- rnorm(10000, mean = 1000, sd = 1)


#散布図を作成し平均を描画
plot(norm1,norm2)
abline(h = mean(norm2),v = mean(norm1),col = "red")

 

まず分散1で適当な正規分布を2つ作成。

 散布図はこうなる。赤線はそれぞれの平均。

f:id:kmooog:20160113105141p:plain

 

ここで相関係数は

 

"散布図に平均点で縦横に軸を描き込んで、プロットされている各人の点からそれぞれの軸に下ろした垂線と軸とで囲まれた長方形の(正負つきの)面積の平均"

相関係数のわかりやすい解釈(続き) - たけみたの脱社会学日記

 

なので適当に外れ値選んで長方形を書く。

 

# 外れ値において長方形作成
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2),mean(norm2),norm2[which.max(norm1)],norm2[which.max(norm1)]), col="gray")

 

f:id:kmooog:20160113105350p:plain

 

これらの分布に相関はないので共分散は0.

 

次にnorm2を(norm1+norm2)に変更し同様の解析。

 

# norm1,norm2+norm1で同様の処理
plot(norm1,(norm2+norm1))
abline(h = mean(norm2+norm1),v = mean(norm1),col="red")
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2+norm1),mean(norm2+norm1),norm2[which.max(norm1)]+norm1[which.max(norm1)],norm2[which.max(norm1)]+norm1[which.max(norm1)]), col="gray")

 

norm1が足された分、すべての点はY軸方向正の方向へ移動し、相関が出てくる。

 

f:id:kmooog:20160113105728p:plain

 

ここで今回の図と最初の図の長方形の面積の差を考える。

最初の図の長方形は

(norm1-norm1の平均)*(norm2-norm2の平均)

二番目の図の長方形は

(norm1-norm1の平均)*{(norm1+norm2)-(norm1の平均+norm2の平均)}

この差分を取るとnorm1^2-norm1の平均^2(=norm1の分散)となる。

 

よって今回の共分散はnorm1の分散と等しく1。

 

ここまではOK。次はnorm1と(norm1+norm2)/√2でやってみる

f:id:kmooog:20160113111535p:plain

おっとそうか図が変わらないのか。元々norm1とnorm2に相関無いから共分散に影響を与えるのはこの場合もnorm1の分散で、面積としてはy軸が1/√2になってるので今回の分散はnorm1の分散/√2 = 1/√2。

 

なんかグラフィカルに理解しないと頭に入ってこない。

確率変数の計算得意になりたい。

 

#適当な正規分布作成
norm1 <- rnorm(10000, mean = 100, sd = 1) 
norm2 <- rnorm(10000, mean = 1000, sd = 1)


#散布図を作成し平均を描画
plot(norm1,norm2)
abline(h = mean(norm2),v = mean(norm1),col = "red")
which.max(norm1)

# 外れ値において長方形作成
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2),mean(norm2),norm2[which.max(norm1)],norm2[which.max(norm1)]), col="gray")


var(norm1,norm2)

# norm1,norm2+norm1で同様の処理
plot(norm1,(norm2+norm1))
abline(h = mean(norm2+norm1),v = mean(norm1),col="red")
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2+norm1),mean(norm2+norm1),norm2[which.max(norm1)]+norm1[which.max(norm1)],norm2[which.max(norm1)]+norm1[which.max(norm1)]), col="gray")



#norm1が足された分、すべての点はY軸方向正の方向へ移動し、相関が出てくる

var(norm1,(norm1+norm2))
sum((norm1-mean(norm1))*( (norm1+norm2)-mean( (norm1+norm2))))/10000



plot(norm1,(norm2+norm1)/1.4142)
abline(h = mean(norm2+norm1)/1.4142,v = mean(norm1),col = "red")
sum((norm1-mean(norm1))*( (norm1+norm2)/1.4142 -mean( (norm1+norm2)/1.4142)))/10000
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean((norm1+norm2)/1.4142),mean(norm1+norm2)/1.4142,(norm2[which.max(norm1)]+norm1[which.max(norm1)])/1.4142,(norm2[which.max(norm1)]+norm1[which.max(norm1)])/1.4142), col="gray")



今日一番の気付き

今日まではpythonスクリプトをqsubで投げる時わざわざ.shファイルを作るようなつまらない人生を送ってきましたが、以下のコマンドで済むことを知ってQOL(クオリティ・オブ・ライフ)が一気に上がりました。

 

qsub -cwd -S $(which python) hogehoge.py

 

 

最近見た色々

 

個人的な評価を1.0~10.0で

 

映画

ゾンビ・クエスト 7.5

オランダのゾンビ映画。ショーンオブザデッド、ゾンビランド系の非低予算エンタメゾンビ映画としてはなかなかの出来。でもちょっと演出の振り切れてなさみたいなのは若干あったかな。

 

グランド・ブタペスト・ホテル  9.0

すごい色が綺麗な映画だった。基本コメディなんだけどヒューマン要素も強くてなんか得した気分になれる。雰囲気としてはアメリの監督の「ミックマック」みたいな感じでしょうか。ストーリーじゃなく、絵面が残る映画は好きだな。良い。

 

なんか続きモノのDVDばっか見てて映画あんまり見てないな

 

DVD

 

伊集院光のばらえてぃ 6.5

最近ずっと見てる。売れてない頃のばいきんぐ小峠がいい感じ。ラジオの感じの伊集院は映像で見ると若干違和感あるが最近慣れてきた。

 

アンタッチャブルのできませんはいいません 7.0

偶然これにも伊集院出てた。アンタッチャブルのシカゴマンゴは中学時代、はじめてはまった深夜ラジオかもしれない。コンビ結成時や下積み時代のエピソードとか、今見るとなかなか感慨深いものがある。

 

MAN vs WILD 8.0

これもマイブーム。元特殊部隊で最年少エベレスト登山に成功したイギリス人のおっちゃんがジャングルで生き延びるためにひたすら蜘蛛やらトカゲやらを捕獲して喰うシリーズ。本当に不味そうに虫を喰うのがいい。

 

小説

 

銀河ヒッチハイク・ガイド 9.5

なんで今までこれ読んでなかったんだろうというぐらい好みのSFだった。久々に時間を忘れて本を読んだ。

 

漫画

 

デッドデッドデーモンズデデデデデストラクション 8.0

浅野いにお食わず嫌いしてたけど面白いじゃないか!

 

夏目友人帳 9.0

言わずと知れた名作。アニメは見てたけど漫画もいいね。全話で泣いてる気がする。

blasrのインストールに再挑戦

環境は遺伝研スパコン

 

LPM: Local Package Managerにパッケージが上がっていたので試してみた。

hdf5は一発で通るのにblasrは通らない。

 

---

iblasr/BlasrHeaders.h:22:23: 致命的エラー: libconfig.h: そのようなファイルやディレクトリはありません

---

致命的エラー。そうですか。

 

仕方ないのでblasrに必要なpbcoreを自力で入れようと試みる。

 

cd pbcore

python setup.py install

 

---

The installation directory you specified (via --install-dir, --prefix, or

the distutils default setting) was:

    /usr/local/lib/python2.7/site-packages/

---

 

デフォルトではpythonsite-packagesに入れることになってるのか?やり直し。

 

python setup.py install  --prefix /home/***

 

---

running install

Checking .pth file support in /home/***/lib/python2.7/site-packages/

error: can't create or remove files in install directory

The following error occurred while trying to add or remove files in the

installation directory:

    [Errno 2] No such file or directory: '/home/***/lib/python2.7/site-packages/test-easy-install-35953.pth'

The installation directory you specified (via --install-dir, --prefix, orthe distutils default setting) was:

    /home/***/lib/python2.7/site-packages/

This directory does not currently exist.  Please create it and try again, or choose a different installation directory (using the -d or --install-dir option).

 

---

なんかsite-packagesに入れたがってるっぽい。

ローカルにpython2.7を落としてsite-packageを指定してやり直し。

 

export PYTHONPATH=/home/***/lib/python2.7/site-packages/

python setup.py install --prefix /home/***

 

通りそうな雰囲気を出しつつも、

ImportError: No module named Cython.Build

で止まる。えぇー

 

 

とりあえずお手上げ。

盗まれた自転車が戻ってきた話

一昨年の冬、自転車が盗まれた。

夕方にラボから戻って、家でちょっと準備してから夕飯の買い物に行こう、というその10数分の間鍵をかけずアパートの駐輪場に置いておいたところを盗まれたのだ。

ママチャリではあったが買ったばっかりだったのでそれなりに悲しく、また今住んでいる地域の治安の悪さを実感して嫌な気持ちになる事件だった。

一応自転車屋さんの盗難保険に入っていたので、警察の人に来てもらい盗難届けを出し、保険でかなり割り引かれた状態で新しい自転車を手に入れられたのでそれっきり忘れていた。

 

しかし3日ほど前、突然電話がかかる。

「もしもし、柏警察署ですがこちら熊谷さんのお電話でよろしかったでしょうか?」

最近警察に迷惑をかけるようなことをした覚えもないし、ここ数年はややこしい人たちとの付き合いもない。即座に自転車が見つかったことだろうと察する。話を聞くと、ごく家の近所で盗まれた自転車が見つかったらしい。

東京に住んでいた頃は、自転車で走っているだけでしょっちゅう呼び止められて登録番号を確認されたものだが、柏に移ってからそのようなことが一切ないので、盗難された自転車が見つかることはまずないだろうと思っていた。柏の外れでも、社会秩序を保とうと努力している人はいるようだ。

「一応盗難事件の証拠品という扱いになりますので、本日現場に来れない場合は柏警察署まで来ていただくことになりますがよろしいでしょうか?」

よろしいも何もそうするしかないんでしょうそうしますよ、と昔さんざん職務質問を食らった経験から来る警察への卑屈さを滲ませながら、後日取りに行く約束をする(その日は学部時代にお世話になった先輩の結婚式の二次会があり、久々に高田馬場の汚い空気を吸いに行くところであった)。

 

「車ないので取りに行くとしたら乗って帰ることになると思うんですが、自転車は使える状態でしょうか?」

「見たところ特に問題ないですよ!全然乗って帰れると思います!」

 

そんな顛末でわざわざ柏警察署まで出向き、30分以上かけて乗って帰ってきた自転車がこちら。

 

f:id:kmooog:20150930152852j:plain

f:id:kmooog:20150930152906j:plain

 

まごうことなき粗大ゴミである。

もしも僕がピクサーアニメの監督で、自転車を擬人化した作品を作っていて、レイプされた後のママチャリを表現しなければいけないなら、この自転車は非常にいいモデルになるのではないか。

バイク風になるように、わざわざブレーキの位置を変えているあたり盗んだ犯人の人柄が透けて見えて好ましい。呼び鈴も外されてるけど、お巡りさんに問題なく乗れると言われたので、道路交通法や条例などの面でもセーフという判断でいいんですね?お巡りさん?これで注意されたらあなたの名前出しますよ?

 

とりあえず、ラボ内の移動用などに自転車欲しかったのでちょっと整備して使う(呼び鈴はつけます)。

最近観た映画

1.0~10.0で書いているのは評価ではなく好みですね

見た順に列挙していきます 

なんかエドガー・ライトばっかり見てた気がする

あと前回書いてからそんなに経ってないのにもう何見たかかなり忘れている

 

DVDなど

5.0 ワールズ・エンド 酔っ払いが世界を救う

「ショーンオブザデッド」「ホットファズ」を加えた三部作の中では一番微妙。主人公に魅力がない。ラストも良くない。酔っ払いvsエイリアンという構図はとても好きなだけに残念。「グラバーズ」という映画が割とコンセプトが似ていて面白そうなのでそっちに期待しよう。

 

8.5 スコット・ピルグリム VS. 邪悪な元カレ軍団

原作のアメコミが好き、かつ監督がエドガー・ライトだったので期待して見たんだけど期待を上回る出来。ビデオゲームをモチーフとした演出が非常にハマっている。主人公、完全にクズなのに憎めないのはなんでなんだ。

4.0 アタック・ザ・ブロック

これもエドガー・ライトが関わってるから見たんだけどエドガー・ライト感はニック・フロストの出演以外に皆無だった。パニック映画としての出来はかなりいいのだけど、この手の映画は借りてまで見るもんじゃないというのが僕の思想なので定評化。なんにもやることない夏の休日の昼間に、たまたまテレ東でやってるのを冷やし中華食いながらだらだら見る、そんな出会い方をしたかった。

 

4.0 アバウトシュミット

みんな大好きジャック・ニコルソン。本作では抑圧された気の弱い平凡なサラリーマンの老後という、およそらしくない役を演じていて、最初は(先入観もあり)ちょっと違和感あるんだけどだんだんそうとしか見えなくなってくる。役者の演技についての感想とか、映画見てほとんど抱くことないけどニコルソン御大は別格だなー

 

5.5 モンスターズ・インク

はじめて見たんだけどピクサー映画の中でとりたてて名作という感じはしなかったなぁ。吹き替え版のウーチャカの声優は相当いい。

blastの出力からヒットした配列を抜き出す

blastをかけて、ヒットした部分の配列をfastaで出力するコード(biopythonとかまともに使えれば絶対こんなのあるんだろうな・・・)

makeblastdbは事前にやっておく

ゲノムからあるタンパク質周辺のシンテニーを知りたい場合などは、周辺領域の長さを指定することでごっそり抜き出したり

 >|python|
#!/usr/bin/env python
#coding: UTF-8

#よく使うモジュールは何も考えずに書く

import sys
import re
import os
import glob
import copy
import commands
from bisect import bisect
from Bio import SeqIO
import urllib2
import datetime


"""
使い方
blastEX blasyn query.fasta db.fasta num(default:0)

で走らせると以下のコマンドを実行し、ヒットした配列を取得

blastn -query query -db db.fasta -outfmt 6

numは取得する周辺配列の長さ。
num = 100 にすればヒットした配列の前後100aaを取得。
前後が100aa以下の場合は取れるだけの長さを取る。

"""

now = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")

argvs = sys.argv # コマンドライン引数を格納したリストの取得
argc = len(argvs) # 引数の個数

blast = argvs[1]

#num変数の設定

if argc > 4:
num = int(argvs[4])
else:
num = 0

#blastの実行

os.system("%s -query %s -db %s -outfmt 6 > BlastResult_%s"%(blast,argvs[2],argvs[3],now))
outputFasta = open("Output_%s.fasta"%now,"w")

#配列の取得

for line in open("BlastResult_%s"%now) :
#print line
if len(line) > 2:
ID = line.split("\t")[1]
seqstart = min(int(line.split("\t")[8]),int(line.split("\t")[9]))
seqend = max(int(line.split("\t")[8]),int(line.split("\t")[9]))
index = 0
handle = open(argvs[3], "rU")
for record in SeqIO.parse(handle, "fasta"):
if record.id == ID:
if seqstart < num and seqend > len(record.seq)-num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq ) +"\n" )
elif seqstart < num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[0:seqend + num] ) +"\n" )
elif seqend > len(record.seq)-num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[seqstart-num:] ) +"\n" )
else:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[seqstart-num:seqend + num]) +"\n" )
handle.close()
outputFasta.close() ||<