怪しい研究成果の発表について（その3）

いうまでもないけど，こういうこともある

 

ビッグデータだ!サンプルサイズは 100,000,000 だ。

変数だって同じくらいあるんだぜ!!

おそれいったか。　（はい，おそれいりました。。震え声）

おそれながら R で例示するのには，100,000,000 は大きすぎるので，1000 で勘弁してください（オロオロ）

> n = 1000 # サンプルサイズ 1000 のデータを
> p = n-1 # 999 個の説明変数で予測する
> y = rnorm(n) # 被説明変数
> x = matrix(rnorm(n*p), n) # 説明変数
> a = lm(y ~ x) # 重回帰分析します
> all.equal(unname(predict(a)), y) # ほらね？予測値は実測値と全く同じ，つまり，100%正確に予測できるんですよ（オッホン）
[1] TRUE

アふぉ

怪しい研究成果の発表について（その2）

わざわざ言う必要もないのだけど，件の研究発表は，以下のようなくだらない解析結果だったのですよ。

> set.seed(1)
> d = data.frame(y = rnorm(20), x = letters[1:20])
> d
             y x
1  -0.62645381 a
2   0.18364332 b
    :
20  0.59390132 t

つまり，データの発生源を独立変数（説明変数）にしてしまった。

分析者は，1個の独立変数を使ったつもりだけど，実際は19個のダミー変数を使ってしまった。20 個のデータポイントを予測するのに，19個の変数を使うと，100%の予測が出来てしまう。それ以外の独立変数を使おうが使うまいが，100%に変わりはない。

> ans = lm(y ~ x, d)
> summary(ans) # 解は求まるが，何の意味もない。独立変数と従属変数は 1:1 の関係値にあるのだから，独立変数がある値を取るときは，対応する従属変数を予測値として返すと言うだけの話。

Call:
lm(formula = y ~ x, data = d)

Residuals:
ALL 20 residuals are 0: no residual degrees of freedom!

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.626454         NA      NA       NA
xb           0.810097         NA      NA       NA
xc          -0.209175         NA      NA       NA
    :
xt           1.220355         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1,    Adjusted R-squared:    NaN
F-statistic:   NaN on 19 and 0 DF,  p-value: NA

> cbind(d$y, predict(ans))
          [,1]        [,2]
1  -0.62645381 -0.62645381
2   0.18364332  0.18364332
3  -0.83562861 -0.83562861
    :
20  0.59390132  0.59390132

> all.equal(d$y,unname(predict(ans))) # 予測値は実測値と完全に一致する!!!
[1] TRUE

＝＝＝＝＝＝＝＝

R など，Factor（カテゴリー） 変数を自動的に処理してくれる統計処理システムだからこそ起こる悲劇。

もし，data.frame の x が自動的に Factor にされないとすると，d$x は

> d$x = as.integer(d$x)

として扱われたのと同じになる。

これだと，

> ans2 = lm(y ~ x, d)
> summary(ans2)

Call:
lm(formula = y ~ x, data = d)

Residuals:
    Min      1Q  Median      3Q     Max
-2.4808 -0.5168  0.1875  0.4833  1.5450

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03609    0.43158  -0.084    0.934
x            0.02158    0.03603   0.599    0.557

Residual standard error: 0.9291 on 18 degrees of freedom
Multiple R-squared:  0.01955,    Adjusted R-squared:  -0.03492
F-statistic: 0.3589 on 1 and 18 DF,  p-value: 0.5566

となるので，逆に疑わしい結果が得られる。

つまり，本来ダミー変数 a，b，c，... であるべきものが 1, 2, 3,  ... として使われるのだ。

この違いを知らないと，地獄に落ちる。

おお，神よ!!（^_^）

怪しい研究成果の発表について

怪しいんですよと解説してくれているサイトでの表示が，奥歯にものが挟まったような気がするので...

まず，データをダウンロードして実行して見ると，件のサイトには明示されていない「警告メッセージ」が現れる。

>  d.glm<-glm(Result~.,d,family=binomial)
 警告メッセージ:
1:  glm.fit: アルゴリズムは収束しませんでした  
2:  glm.fit: 数値的に 0 か 1 である確率が生じました  

つまり，「以下に示される結果は『信頼できませんよ』ということですかね

それを無視して，分析を進めると

> table(d$Result,round(predict(d.glm,newdata=d[,-3],type='response'),0))
   
      0   1
  0 251   0
  1   0 240
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading

またしても，警告メッセージが出ていますけど。（表示される結果は，件のサイトと同じ）

以下，件のサイとで行われたとおりの分析を続けると，出てくる結果はそのサイトに示されているのと全く同じ。

違う所（というか，件のサイとの結果の表示では省略されているだけかも知れないのだけれど），「警告メッセージ」というのが出てくるのよ。

こんなのが，毎回出てくると「心穏やかにならなくなるよね」

> d.glm.true<-glm(Result~.,d[,-c(1,2)],family=binomial)
> table(d$Result,round(predict(d.glm.true,newdata=d[,-c(1,2,3)],type='response'),0))
   
      0   1
  0 238  13
  1  14 226
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading
> library(e1071)
> d.svm<-svm(Result~.,d,kernel='linear')
> table(d$Result,predict(d.svm,newdata=d[,-3]))
   
      0   1
  0 251   0
  1   1 239
> d.svm.true<-svm(Result~.,d[,-c(1,2)],kernel='linear')
> table(d$Result,predict(d.svm.true,newdata=d[,-c(1,2,3)]))
   
      0   1
  0 238  13
  1  15 225
> nrow(d[,1:2])
[1] 491
> nrow(unique(d[,1:2]))
[1] 485
> d.glm.name<-glm(Result~.,d[,1:3],family=binomial)
 警告メッセージ:
1:  glm.fit: アルゴリズムは収束しませんでした  
2:  glm.fit: 数値的に 0 か 1 である確率が生じました  
> table(d$Result,round(predict(d.glm.name,newdata=d[,1:2],type='response'),0))
   
      0   1
  0 247   4
  1   2 238
 警告メッセージ:
 predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  で:
  prediction from a rank-deficient fit may be misleading
> summary(d.glm)

Call:
glm(formula = Result ~ ., family = binomial, data = d)

Deviance Residuals:
       Min          1Q      Median          3Q         Max  
-7.228e-06  -2.409e-06  -2.110e-08   2.409e-06   8.367e-06  

Coefficients: (47 not defined because of singularities)
                                Estimate Std. Error z value Pr(>|z|)
(Intercept)                    1.184e+02  4.986e+06       0        1
Player1A.Kuznetsov            -3.999e+01  1.553e+06       0        1
Player1A.Man0rino             -2.280e+01  1.667e+06       0        1
Player1A.Ramos                -1.334e+01  2.653e+06       0        1
Player1A.Seppi                 2.574e+01  1.920e+06       0        1
Player1A.Ungur                 6.009e+01  1.841e+06       0        1
Player1Adrian Man0rino        -2.082e+01  2.717e+06       0        1
　: あまりにも多いので，途中省略

BPC.2                         -6.438e+00  9.798e+04       0        1
BPW.2                          1.426e-01  4.375e+04       0        1
NPA.2                         -1.442e-01  6.894e+04       0        1
NPW.2                          4.440e-01  7.039e+04       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6.8042e+02  on 490  degrees of freedom
Residual deviance: 3.9260e-09  on  74  degrees of freedom
AIC: 834

Number of Fisher Scoring iterations: 25

この結果見ただけで，普通の人は，「げんなりする」はず。

原論文の解析者は，この結果を，見ていないのかな（みていないんだろうな）。

ビッグダーどころか，統計学，統計データ，多変量解析について，なんの知識もない人なんだなあというのが，正直な感想（卒論レベルの評価でも，不合格ですよ）。筆者は工学分野が専門なようで，当該のような医学・薬学分野の常識が欠除していたんだろうなあと思います。

このような人を特任研究員？として任用した，京都大学の然るべき部署の責任が問われるべきでしょうね。

一度は公表して，すぐに取り下げたということであるが，責任ある京都大学の然るべき部署として，なぜ取り下げたか，どういう根拠に基づいて取り下げるを得なかったのかをしっかりと説明しないと，ダメですよ～～～～～。って。嘆かわしい。

攪乱順列

https://codeiq.jp/q/2636

銀行のATMなどで暗証番号を入力するときに使用するテンキー。
覗き見られたときのセキュリティを意識して、配置がランダムに変わる機種もあります。
今回はすべてのキーが元の配置とは異なる位置に配置することになりました。
そこで、このような配置が何通りあるかを求めるように依頼されました。

ビット演算だとかメモ化だとか再帰プログラムとかいっている。

「メモ化の練習」と言っているのだからまあよいが，そんなもん，公式があるんだから，瞬殺。

http://mathtrain.jp/montmort

これに限らないが，なんでもプログラムでやろうとするのは止めといた方がよい。

> func = function(n) factorial(n) * sum((-1)^(2:n) / factorial(2:n)) # 一行野郎

> system.time(print(func(10)))
[1] 1334961
   ユーザ   システム       経過  
         0          0          0
> system.time(print(func(12)))
[1] 176214841
   ユーザ   システム       経過  
         0          0          0

n = 18 までは正しい答えが得られる

gawk で --bignumオプションを使って最適化すると，無制限

$ cat b.awk
{
  n = $1
  sign = s = n % 2 == 0 ? 1 : -1
  p = 1
  for (i = n; i >= 3; i--) {
    p = p*i
    sign = -sign
    s = s + sign * p
  }
  printf "%i\n", s
}

$ echo 10 | gawk --bignum --file=b.awk
1334961

$ echo 12 | gawk --bignum --file=b.awk
176214841

$ echo 100 | gawk --bignum --file=b.awk
34332795984163804765195977526776142032365783805375784983543400282685180793327632432791396429850988990237345920155783984828001486412574060553756854137069878601

$ echo 101 | gawk --bignum --file=b.awk
3467612394400544281284793730204390345268944164342954283337883428551203260126090875711931039414949888013971937935734182467628150127669980115929442267844057738700

2次元の数列（発展系）

以下のような 12×12 の行列がある。

          1          1          2          4          7         13         24         44         81        149        274        504
          2          3          7         16         33         69        142        288        580       1159       2301       4544
          3          7         19         49        115        265        595       1307       2828       6038      12748      26662
          5         15         46        131        340        851       2059       4845      11163      25264      56321     123954
          8         30        103        321        909       2449       6333      15843      38616      92094     215622     496948
         13         58        220        743       2270       6533      17938      47429     121679     304404     745455    1792440
         21        109        453       1647       5388      16470      47776     132906     357447     934627    2386057    5967519
         34        201        908       3533      12300      39744     121291     353670     993831    2707823    7186836   18648449
         55        365       1781       7381      27215      92591     296254     902636    2642759    7484099   20602387   55345213
         89        655       3433      15091      58694     209553     700883    2225436    6772462   19890703   56677824  157334651
        144       1164       6522      30302     123897     462865    1614201    5329035   16821322   51139360  150569928  431210474
        233       2052      12240      59918     256801    1001377    3633180   12445829   40674170  127783242  388150993 1145153530

同じ規則で作られる行列の 100行100列目の数値はいくつか？

GAWK で mpfr を使いたいんだが...（その2）

古い Macintosh （Mac OS X 10.7.5）にインストールしてみたら，何の問題もなくインストールできた。

新しい Macintosh の方のディレクトリー構造がおかしくなっているのかな？

$ gawk --version
GNU Awk 4.1.3, API: 1.1 (GNU MPFR 3.1.3, GNU MP 6.1.0)

ちゃんと計算できている

$ gawk --bignum "BEGIN{print 5999856**3+99992800**3-100000000**3}"
-2985984

$ gawk --bignum "BEGIN{a=1;for (i = 1 ; i <= 60; i++) a = a*i; print a}"
8320987112741390144276341183223364380754172606361245952449277696409600000000000000

$ gawk --bignum "BEGIN{print 2**300}"
2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397376

$ gawk --bignum "BEGIN{n=400;a=b=1;for (i = 3; i <= n; i++) { c = a+b; a=b; b=c};print c}"
176023680645013966468226945392411250770384383304492191886725992896575345044216019675

GAWK で mpfr を使いたいんだが...

gmp と mpfr をインストールした状態で gawk をインストールすると，それらが使えるというんだが，できないんだなあこれが。

だれか，教えて。

＝＝＝

失礼。詳細を書かなかった。

Mac OS X El Capitan 10.11.2

gawk, gmp, mpfr は最新版をソースから make

gawk-4.1.3
gmp-6.1.0
mpfr-3.1.3

どれも，make check では問題は見つからないが，

$ gawk -M
gawk: warning: -M ignored: MPFR/GMP support not compiled in

となるのだね。

「電卓」の仕様（その3）

「ほどほどの大きさの整数を扱える，足し算，引き算，掛算，累乗ができる。ただし，挿入演算子を使うのではない」という仕様の電卓プログラム

digits = 100
carry = function(res) {
    len = length(res)
    cry = 0
    for (i in 1:len) {
        t = res[i] + cry
        res[i] =t %% 10
        cry = t %/% 10
    }
    if (cry) res[len+1] = cry
    res
}

add = function(a, b) {
    if (length(a) == 1 && class(a) != "LongInt") a = as.LongInt(a)
    if (length(b) == 1 && class(b) != "LongInt") b = as.LongInt(b)
    if (a[digits] == 9) a = -compliment(a)
    if (b[digits] == 9) b = -compliment(b)
    res = carry(a+b)[1:digits]
    class(res) = "LongInt"
    res
}

sub = function(a, b) add(a, -b)

compliment = function(a) {
    add(9-a, c(1, rep(0, digits-1)))[1:digits]
}

operand = function(a) {
    if (length(a) == 1 && class(a) != "LongInt") a = as.LongInt(a)
    sign = 1
    if (a[digits] == 9) {
        a = compliment(a)
        sign = -1
    }
    return(list(sign= sign, operand=cutZero(a)))
}

mult = function(a, b) {
    a0 = operand(a)
    a = a0$operand
    b0 = operand(b)
    b = b0$operand
    len.a = length(a)
    len.b = length(b)
    res = integer(len.a+len.b)
    for (j in 1:len.b) {
        for (i in 1:len.a) {
            res[i+j-1] = res[i+j-1] + a[i]*b[j]
        }
    }
    res = fillZero(carry(res))
    if (length(res) > digits) {
        stop("Overflow")
    } else if (a0$sign * b0$sign == -1) {
        res = compliment(res)
    }
    class(res) = "LongInt"
    res
}

cutZero = function(a) {
    a = rev(a)
    res = rev(a[cumsum(a) > 0])
    res
}

fillZero = function(a) {
    res = a
    len = length(a)
    if (digits > len) res = c(res, rep(0, digits-len))
    res
}

as.LongInt = function(a) {
    res = rev(unlist(strsplit(as.character(a), "")))
    len = length(res)
    if (res[len] == "-") {
        res = compliment(fillZero(as.integer(res[-len])))
    } else {
        res = as.integer(res)
    }
    res = fillZero(res)
    class(res) = c("LongInt", "integer")
    res
}

print.LongInt = function(a) {
    if (all(a == 0)) {
        a = 0
    } else {
        Sign = ""
        if (a[digits] == 9) {
            a = compliment(a)
            Sign = "-"
        }
        a = paste0(Sign, paste(rev(cutZero(a)), collapse=""))
    }
    cat(a, "\n", sep="")
    return(a)
}

power = function(a, n) {
    if (length(a) == 1 && class(a) != "LongInt") a = as.LongInt(a)
    result = as.LongInt(1)
    for (i in 1:n) {
        result = mult(result, a)
    }
    result
}
square = function(a) power(a, 2)
cubic = function(a) power(a, 3)

###

> add(-13, 3)
-10
> mult(123, 5)
615
> sub(add(square(3), square(4)), square(5))
0
>
> a = as.LongInt(5999856)
> b = as.LongInt("99992800") # 大きな数値は文字型で与えると安心
> c = as.LongInt("100000000")
> sub(add(cubic(a), cubic(b)), cubic(c))
-2985984
>
> # 階乗
> a = as.LongInt(1) # これは必須
> for (i in 1:60) a = mult(a,i)
> a
8320987112741390144276341183223364380754172606361245952449277696409600000000000000
> library(gmp)
> factorialZ(60)
Big Integer ('bigz') :
[1] 8320987112741390144276341183223364380754172606361245952449277696409600000000000000
>
> # 2 のべき乗
> power(2, 300)
2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397376
> as.bigz(2)^300
Big Integer ('bigz') :
[1] 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397376
>
> # フィボナッチ数列
> n = 400
> a = as.LongInt(1)
> b = as.LongInt(1)
> for (i in 3:n) {
+     c = add(a, b)
+     a = b
+     b = c
+ }
> c
176023680645013966468226945392411250770384383304492191886725992896575345044216019675
> fibnum(n)
Big Integer ('bigz') :
[1] 176023680645013966468226945392411250770384383304492191886725992896575345044216019675

「電卓」の仕様（その2）

Java に BigDecimal というのがあるのを知っていたので，どうせ BigInteger というのもあるだろうと，門前の小僧が習わぬ経を詠んでみたという話

ideon.com でやってみた。試行錯誤の末。

整数方程式

締め切りが 2016/1/14 10:00 AM なので，その 1 分後に公開されるように予約投稿する

2つの自然数の組 (a, b) が与えられたとき、自然数 x, y に関する次の方程式を考えます。
　　　　x2 ＋ a2 ＝ y2 + b2　…　(※)
例えば、(a, b) = (3, 10) のとき、方程式(※)の解は (x, y) = (10, 3), (46, 45) の 2 組です。
自然数の組 (a, b) に対し、方程式(※)の全ての解の x + y の和を F(a, b) と定義します。
例えば F(3, 10) = 10 + 3 + 46 + 45 = 104 です。
同様に、F(10, 50) = 3500、F(20, 100) = 15022 となることが確かめられます。
標準入力から、半角空白区切りで 2 つの自然数 a, b（1 ≦ a < b ≦ 105）が与えられます。
標準出力に F(a, b) の値を出力するプログラムを書いてください。

探索範囲を効率的に定めて，ブルートフォースで片付ける

f = function(a, b) {
    s = 0
    b2a2 = b^2-a^2
    xmy = 2 - (b2a2 %% 2)
    repeat {
        xpy = b2a2 %/% xmy
        if (xmy*xpy == b2a2) {
            x = (xmy+xpy)/2
            y = xpy-x
            if (x == floor(x) && y > 0) {
                s = s+x+y
            }
        }
        if (xmy >= xpy) break
        xmy = xmy+2
    }
    cat(s)
}

f(3, 10) # 104
f(10, 50) # 3500
f(20, 100) # 15022
f(10, 26) # 716
f(11, 389) # 291886
f(123, 456) # 294166
f(35672, 61243) # 5082405064
f(71200, 82321) # 2490724704
f(19126, 98765) # 16043311668

「電卓」の仕様

こういうバグもあった。 フェルマーの最終定理の「反例」（電卓編）

http://blog.unfindable.net/archives/621 

ということであるが，ここでいっている「電卓」は，コンピュータ上に実装されているもの。

「電卓の」仕様はさまざまというか，単純に倍精度実数を使って実装されているならば，まさにバグでもなんでもない。Windows の電卓アプリは素晴らしくて，Mac のはダメだというのは筋違い。

99992800^3 は 10 進数で何桁になるか?

> 3*log10(99992800)
[1] 23.99991

そう。24桁だ。dobule で正確に表現できる整数は たかだか 15,6 桁なので，そんな大きな数は扱いきれない。（言うまでもないことであるが 15,6 桁の範囲内の整数であれば，「コンピュータは2進数で表現するので，誤差がある」なんてことはない。整数は常に「正確に」表される）。有効数字が足りない。15，6桁目以降はゴミだ。

このブログの著者は，

> 9999999999999999999999999^9999999
[1] Inf

という結果が出るのもバグだというのだろうかね。

扱いきれない，大きな2数の和から大きな数を引いて，正確な答えが出るわけがない。だからダメだなんて言ってもお門違いも甚だしい。

正しい答えを得るにはどうするか。大きな数を扱える仕様を策定しないといけない。

簡単な四則演算関数を自作するもよし，すでにあるパッケージ（gmp）を使うも良し。

Rmpfr は，使用するビット数を指定できるので，ある桁数の数字を表現するのに何ビット必要か体感できる。

> library(Rmpfr)
> prec = 80 # double の精度なら　53 だ。
> a = mpfr("5999856", prec)
> b = mpfr("99992800", prec)
> c = mpfr("100000000", prec)
> print(a^3+b^3-c^3)
1 'mpfr' number of precision  80   bits
[1] -2985984