Last active
August 29, 2015 14:25
-
-
Save evanbiederstedt/a763375f068b05167c26 to your computer and use it in GitHub Desktop.
m^t C^-1 m, vary C3 from e-08 to e-12, hold C3
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Same as https://gist.github.com/evanbiederstedt/72b0035ca4b9fce830c3 | |
except vary C3, hold sigma^2 fixed | |
FIRST, NO PEAK CODE, i.e. scale covariance matrix by e+21 | |
CODE | |
"""" | |
# | |
# Hold sigma^2 constant, vary C3 | |
# | |
vary_x_samples125 = np.logspace(-8, -12, num=40) # C3 parameter, vary from e-08 to e-12 | |
sigma125 = 5e-22 # chose this sigma^2 parameter, hold constant | |
Sij = vary_x_samples125[:, None, None] * norm_matrix[1][None, :, :] | |
newSij = (1e21)*Sij # multiply S_ij by 1e12 | |
Nij = sigma125 * id_mat[None, :, :] | |
newNij = (1e21)*Nij | |
# Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l | |
# Sij.shape = (40, 3072, 3072) | |
Cij = newSij + newNij | |
#invCij = np.linalg.inv(Cij) | |
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1] | |
# model_fit_terms = m^T C^-1 m | |
# | |
# model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) | |
# for i in range(invCij.shape[0])]) | |
# | |
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) | |
print "m^T c^-1 m terms are" | |
print model_fit_terms | |
print "******************" | |
print "logdetC terms are" | |
print logdetC[1] | |
print "******************" | |
print "Npix2pi term is" | |
print Npix2pi | |
"""" | |
OUTPUT: | |
m^T c^-1 m terms are | |
[ 17424944.59432721 25431124.82535299 20988175.41728179 18788639.9688897 | |
17380549.27187007 16487860.85815524 15909766.35792065 | |
15497517.16395973 15189182.91137147 14972008.39197645 | |
14793336.73669815 14664017.69568876 14568322.70447167 14492458.7918589 | |
14433833.23682215 14389138.7429928 14352751.6784684 14325011.10904521 | |
14302605.58943325 14285747.28930357 14272097.84990945 | |
14261432.05420008 14253089.98221404 14246570.26856775 | |
14241465.54843466 14237157.87713204 14234018.59879179 | |
14231438.17048982 14229450.48704928 14227802.58328702 | |
14226586.63568589 14225606.63455151 14224772.51924103 | |
14224178.22232257 14223689.86871848 14223295.19405429 | |
14222996.47245571 14222746.69923505 14222557.92339739 | |
14222409.72199761] | |
****************** | |
logdetC terms are | |
[-1910.79269502 -1899.49346143 -1894.31400072 -1893.9865907 -1890.0321336 | |
-1889.55785882 -1889.79498435 -1890.52858638 -1891.04930553 -1892.56437866 | |
-1893.64297778 -1895.06109716 -1896.61902812 -1898.36601353 -1899.78459156 | |
-1901.39518245 -1903.05163531 -1904.67733691 -1906.32176722 -1907.95386497 | |
-1909.59522293 -1911.21625436 -1912.86998408 -1914.50988852 -1916.15690315 | |
-1917.81120861 -1919.462296 -1921.10832586 -1922.76238417 -1924.41349353 | |
-1926.06744684 -1927.71856858 -1929.37104233 -1931.02466563 -1932.67642731 | |
-1934.32805631 -1935.98147704 -1937.63429343 -1939.28755537 -1940.94054106] | |
****************** | |
Npix2pi term is | |
19301.9452637 | |
************************************** | |
************************************* | |
SECOND, THERE IS A PEAK CODE, i.e. scale covariance matrix by e+22 | |
CODE | |
"""" | |
# | |
# Hold sigma^2 constant, vary C3 | |
# | |
vary_x_samples123 = np.logspace(-8, -12, num=40) # C3 parameter, vary from e-08 to e-12 | |
sigma123 = 5e-22 # chose this sigma^2 parameter, hold constant | |
# | |
# Plot noise sigma from e-22 to e-24 | |
# | |
Sij = vary_x_samples123 [:, None, None] * norm_matrix[1][None, :, :] | |
newSij = (1e22)*Sij # multiply S_ij by 1e12 | |
Nij = sigma123 * id_mat[None, :, :] | |
newNij = (1e22)*Nij | |
# Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l | |
# Sij.shape = (40, 3072, 3072) | |
Cij = newSij + newNij | |
#invCij = np.linalg.inv(Cij) | |
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1] | |
# model_fit_terms = m^T C^-1 m | |
# | |
# model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) | |
# for i in range(invCij.shape[0])]) | |
# | |
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) | |
print "m^T c^-1 m terms are" | |
print model_fit_terms | |
print "******************" | |
print "logdetC terms are" | |
print logdetC[1] | |
print "******************" | |
print "Npix2pi term is" | |
print Npix2pi | |
"""" | |
OUTPUT | |
m^T c^-1 m terms are | |
[ 2931868.99027531 2561423.09736563 2092876.99077943 1876397.00384942 | |
1740382.76555146 1650113.62054166 1589940.751404 1549446.41704495 | |
1519557.92056381 1496555.37291566 1479682.09032262 1466653.76019673 | |
1456837.51067493 1449220.46772181 1443528.73129329 1438909.39635894 | |
1435211.46326468 1432468.39780895 1430277.65089106 1428567.73979521 | |
1427238.59379346 1426165.74999515 1425328.54426065 1424671.52510419 | |
1424138.61076047 1423723.85218339 1423394.01406752 1423137.00097182 | |
1422944.56767479 1422779.35144192 1422656.21763844 1422559.03605542 | |
1422480.3001671 1422414.48572627 1422367.35435307 1422328.77614777 | |
1422298.52139999 1422275.64094104 1422255.88586292 1422241.38234343] | |
****************** | |
logdetC terms are | |
[ 5162.11916985 5174.89668376 5179.54501215 5178.57071109 5182.81051551 | |
5184.46873299 5183.19398907 5183.19703352 5182.37013867 5181.01067721 | |
5179.97913256 5178.52226416 5176.80653401 5175.16623507 5173.76193185 | |
5172.11691893 5170.45643259 5168.87473124 5167.23262542 5165.6055479 | |
5163.95696043 5162.3116691 5160.67301514 5159.02720484 5157.37926926 | |
5155.73010752 5154.07978174 5152.43292945 5150.77874406 5149.12856074 | |
5147.47438028 5145.82296635 5144.17003878 5142.51772009 5140.86502619 | |
5139.21290009 5137.56004832 5135.9069266 5134.25393896 5132.60101615] | |
****************** | |
Npix2pi term is | |
19301.9452637 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
FIRST, FIND NO PEAK IN LF, NOISE PARAMETERS E-22 TO E-24
CODE, C3 = 5e-10:
""""
Hold C3 fixed, then do test again with 10% less C3 value
Use C3=5e-10
Noise is now varied from e-22 to e-24
vary_x_samples123 = 5e-10
sigma123 = np.logspace(-22, -24, num=40)
Sij = vary_x_samples123 * norm_matrix[1][None, :, :]
newSij = (1e22)_Sij # multiply S_ij by 1e12
Nij = sigma123[:, None, None] * id_mat[None, :, :]
newNij = (1e22)_Nij
Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
Sij.shape = (40, 3072, 3072)
Cij = newSij + newNij
invCij = np.linalg.inv(Cij)
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1]
model_fit_terms = m^T C^-1 m
model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) )
for i in range(invCij.shape[0])])
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
print "m^T c^-1 m terms are"
print model_fit_terms
print "_"
print "logdetC terms are"
print logdetC[1]
print "_"
print "Npix2pi term is"
print Npix2pi
""""
OUTPUT
m^T c^-1 m terms are
[ 7.99112774e+06 9.15862622e+06 1.05174464e+07 1.21218540e+07
1.40481199e+07 1.63849383e+07 1.94345938e+07 2.28107290e+07
2.75426116e+07 3.39869501e+07 4.30589812e+07 5.86874289e+07
8.17059720e+07 1.64401397e+08 -3.15809360e+10 -1.55530576e+08
-7.00025652e+07 -3.69097741e+07 -2.01136741e+07 5.33705542e+06
4.54644402e+06 1.50969154e+07 3.64295994e+07 7.59843057e+07
-5.50891593e+07 1.77650620e+07 5.68860609e+07 1.24610626e+08
-1.18748462e+08 1.00196018e+08 3.26506311e+07 -5.96827008e+06
-9.43865379e+07 -2.91770949e+07 1.05304583e+08 -4.74799471e+07
-1.85701131e+07 1.86812899e+07 -4.27639327e+07 -5.13285598e+07]
logdetC terms are
[ 239.71699157 -122.71079135 -486.33852385 -849.25504606
-1212.4921126 -1576.26351061 -1944.85327619 -2305.753459
-2669.13822474 -3030.65420367 -3396.41771088 -3763.82736171
-4131.60575246 -4506.35071982 -4889.11466545 -5248.4037388
-5612.4490359 -5995.67825433 -6385.48032885 -6779.9332622
-7172.39370076 -7581.395549 -8041.60387949 -8517.11067861
-9002.83276846 -9446.49442793 -9839.83647962 -10172.92306214
-10420.76829252 -10623.50283171 -10792.24792309 -10928.10844176
-11027.19559602 -11100.06733525 -11170.65063404 -11215.66540551
-11261.93661065 -11284.43614935 -11306.60681018 -11326.08454794]
Npix2pi term is
19301.9452637
SECOND CODE,
C3 = 4.5e-10
Hold C3 fixed, then do test again with 10% less C3 value
Now use C3=4.5e-10
Noise is now varied from e-22 to e-24
vary_x_samples123 = 4.5e-10
sigma123 = np.logspace(-22, -24, num=40)
Sij = vary_x_samples123 * norm_matrix[1][None, :, :]
newSij = (1e22)_Sij # multiply S_ij by 1e12
Nij = sigma123[:, None, None] * id_mat[None, :, :]
newNij = (1e22)_Nij
Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
Sij.shape = (40, 3072, 3072)
Cij = newSij + newNij
invCij = np.linalg.inv(Cij)
logdetC = np.linalg.slogdet(Cij) # returns sign and determinant; use logdetC[1]
model_fit_terms = m^T C^-1 m
model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) )
for i in range(invCij.shape[0])])
model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
print "m^T c^-1 m terms are"
print model_fit_terms
print "_"
print "logdetC terms are"
print logdetC[1]
print "_"
print "Npix2pi term is"
print Npix2pi
OUTPUT
m^T c^-1 m terms are
[ 7.88230937e+06 9.00401772e+06 1.03091948e+07 1.18496346e+07
1.36733594e+07 1.58717504e+07 1.85223083e+07 2.29277071e+07
2.58774395e+07 3.11601583e+07 3.85718641e+07 4.95209561e+07
6.71862827e+07 9.81421763e+07 2.00515338e+08 1.79149266e+08
-1.57174217e+08 -7.01207170e+07 -3.93179729e+07 -2.01955132e+07
-6.98201421e+06 5.88648498e+06 1.83612716e+07 3.94038299e+07
5.39473080e+07 3.01503865e+07 -3.34339516e+07 -1.31643012e+08
-1.18775863e+08 -8.63387028e+06 2.02255477e+08 -5.89407096e+07
-1.31132724e+08 2.82823527e+07 8.34855351e+06 1.41266857e+08
1.68945314e+08 -8.34350175e+07 -4.03057940e+07 1.36690513e+09]
logdetC terms are
[ 238.87906146 -123.59978549 -486.97635053 -849.54919459
-1212.63010956 -1576.01604717 -1942.33073337 -2309.73321034
-2670.91431763 -3031.09627919 -3396.32588831 -3762.14204865
-4128.07960518 -4502.34341668 -4872.33980536 -5255.02597355
-5607.79232715 -5986.94061374 -6370.62068917 -6753.73665052
-7147.26442618 -7537.878678 -7981.39859695 -8433.45165572
-8907.83987127 -9384.80250567 -9848.92410084 -10252.021386
-10558.32732388 -10804.17249772 -11029.95815536 -11179.59272393
-11324.19494314 -11416.3809462 -11503.90735819 -11568.77486146
-11622.48104156 -11656.23132841 -11687.85082231 -11713.95931238]
Npix2pi term is
19301.9452637