print("Hello world.")

#まだ解けていない問題一覧
#
#問1.サイコロを7回投げて、2以上の目が出る確率。
#こたえ:  5/6
#
#問2.確率変数Uが平均0、分散1の正規分布にしたがうときに、P(0.3+2x+U>0)=0.95となるxを求める。
#こたえ: 0.67
#
#問3.0から1までの一様乱数の分散。

ans = 0.3+2*(-0.67)+rnorm(10000)
length(ans[ans>0]/10000)

ans = 0.3+2*(0.67)+rnorm(10000)
length(ans[ans>0]/10000)

#問1.サイコロを7回投げて、2以上の目が出る確率をシミュレーションで求めよ。
x=runif(10^5, min=1, max=7); #print(x)
print(length(x[x>=2])/10^5)

#問1-2.サイコロを7回投げて、2以上の目が出る確率を求めよ。
print(5/6) # サイコロを何回投げようと2以上の目が出る確率は等しい(ベルヌーイ試行のため)。

#問2.確率変数Uが平均0、分散1の正規分布にしたがうときに、P(0.3+2x+U>0)=0.95となるxを求めよ。

# qnorm(0.95)=1.64より、0.3+2x+U=1.64を解く。
# x=(1.34-U)/2=(1.34-N(0,1))/2=N(-1.34,1)/2=N(-0.67,1); x=-0.67として、

x=0.67; U=rnorm(100000,0,1);
ans=0.3+2*x+U>0; print(length(ans[ans==TRUE])/100000)

#qnorm(0.95)=1.64より、0.3+2x+U=1.64をxについて解いてx=(1.34-U)/2=(1.34-N(0,1))/2=-N(-1.34,1)/2=N(-0.67,1)とする。
#x=-0.67のとき、P(0.3+2x+U>0)=0.95となる確率が最大となる。

#数値計算例:
x=0.67; U=rnorm(100000,0,1);
ans=0.3+2*x+U>0; print(length(ans[ans==TRUE])/100000) # 0.95となる。

x=0.67; U=rnorm(100000,0,1);
ans=0.3+2*x+U>0; print(length(ans[ans==TRUE])/100000) # 0.95となる。

#問3. 0から1までの一様乱数の分散を求めよ。
X=runif(10000,min=0,max=1); var(X) # 0.084
X=runif(1000000,min=0,max=1); var(X) # 0.084
X=runif(10000000,min=0,max=1); var(X) # 0.084 (なぜか。)

print("Sage test")

# Hadamard matrix
m=matrix([[1,1],[1,0]]); print(m)
print(det(m)); print(m^(-1)); print(exp(m).n())

import sage.interfaces.gap
sage.interfaces.gap.gap_cmd = "/usr/local/bin/gap"

print(SymmetricGroup(3));
# R = gap.PolynomialRing('Rationals', 2)  ERROR

a=[0,1,2,3,4,5]; # initial cycle
b=[a[0],a[1],a[2],a[3],a[4],a[5]]
a=[b[0],b[1],b[2],b[3],b[4],b[5]]
print(a); print(b)

# ab = ba

x=[0,1,2,3,4,5]; # initial cycle
a=[x[2],x[1],x[0],x[5],x[4],x[3]]
b=[a[a[2]],a[a[1]],a[a[0]],a[a[5]],a[a[4]],a[a[3]]]
c=[x[a[2]],x[a[1]],x[a[0]],x[a[5]],x[a[4]],x[a[3]]]

print(a); print(b); 
# a = [2, 1, 0, 5, 4, 3]
# b = [2, 1, 0, 5, 4, 3]
# RESULT: a is idempotent.

print(c)
# c = [0, 1, 2, 3, 4, 5]. xa=x

print([a[x[0]],a[x[1]],a[x[2]],a[x[3]],a[x[4]]]) # ax=(2,1,0,5,4)
print([x[a[0]],x[a[1]],x[a[2]],x[a[3]],x[a[4]]]) # xa=(2,1,0,5,4)

# Consideration
# ax = xa holds if x is an initial cycle.

b=[5,4,3,2,1,0] # Define a permutation b.

print([a[b[0]],a[b[1]],a[b[2]],a[b[3]],a[b[4]],a[b[5]]]) # ab = [3, 4, 5, 0, 1, 2]
print([b[a[0]],b[a[1]],b[a[2]],b[a[3]],b[a[4]],a[b[5]]]) # ba = [3, 4, 5, 0, 1, 2]

# ab = ba holds.

#a = [3,4,2,0,1]. How can you deduce the inverse of permutation a such that ax = [0,1,2,3,4]?

# Recall the definition of a permutation cycle notation. a is:

# a = [
#ia = [
#ia = [

x=c(25,10,8,7,41,69,24,16,18,38,9,10)
y=c(15,21,7,6,46,64,22,18,23,32,11,9)
sum((x-y)^2/y)

qchisq(0.975,df=6)

# Generate Normal Gaussian

x = runif(10000,0,100000000); #print(x)
mean(exp(-x^2)) # = 0
sum((exp(-x^2)-1)^2)/1000 # = 1

# Define v as normal density
v = exp(-x^2)

sum(v^2+v^2 < 1) # 1000 ~ 100%
table(v^2 + v^2 < 1)