from sklearn import linear_model
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
from matplotlib import style
import numpy as np
import pandas as pd
from statistics import mean
df = pd.read_csv("data/weight-height.csv")
df
| Gender | Height | Weight | |
|---|---|---|---|
| 0 | Male | 73.847017 | 241.893563 |
| 1 | Male | 68.781904 | 162.310473 |
| 2 | Male | 74.110105 | 212.740856 |
| 3 | Male | 71.730978 | 220.042470 |
| 4 | Male | 69.881796 | 206.349801 |
| ... | ... | ... | ... |
| 9995 | Female | 66.172652 | 136.777454 |
| 9996 | Female | 67.067155 | 170.867906 |
| 9997 | Female | 63.867992 | 128.475319 |
| 9998 | Female | 69.034243 | 163.852461 |
| 9999 | Female | 61.944246 | 113.649103 |
10000 rows × 3 columns
male_df = df[df['Gender'] == 'Male']
male_df
| Gender | Height | Weight | |
|---|---|---|---|
| 0 | Male | 73.847017 | 241.893563 |
| 1 | Male | 68.781904 | 162.310473 |
| 2 | Male | 74.110105 | 212.740856 |
| 3 | Male | 71.730978 | 220.042470 |
| 4 | Male | 69.881796 | 206.349801 |
| ... | ... | ... | ... |
| 4995 | Male | 68.860062 | 177.131052 |
| 4996 | Male | 68.973423 | 159.285228 |
| 4997 | Male | 67.013795 | 199.195400 |
| 4998 | Male | 71.557718 | 185.905909 |
| 4999 | Male | 70.351880 | 198.903012 |
5000 rows × 3 columns
male_df = df[df['Gender'] == 'Male'][:200]
male_df
| Gender | Height | Weight | |
|---|---|---|---|
| 0 | Male | 73.847017 | 241.893563 |
| 1 | Male | 68.781904 | 162.310473 |
| 2 | Male | 74.110105 | 212.740856 |
| 3 | Male | 71.730978 | 220.042470 |
| 4 | Male | 69.881796 | 206.349801 |
| ... | ... | ... | ... |
| 195 | Male | 72.294743 | 208.474866 |
| 196 | Male | 67.253325 | 180.779780 |
| 197 | Male | 75.944460 | 231.924749 |
| 198 | Male | 66.316232 | 170.593858 |
| 199 | Male | 71.156053 | 188.271456 |
200 rows × 3 columns
male_df['Height'] = male_df['Height'].apply(lambda x: x * 2.54)
male_df['Weight'] = male_df["Weight"].apply(lambda x: x * 0.45359237)
male_df
| Gender | Height | Weight | |
|---|---|---|---|
| 0 | Male | 187.571423 | 109.721075 |
| 1 | Male | 174.706036 | 73.622792 |
| 2 | Male | 188.239668 | 96.497629 |
| 3 | Male | 182.196685 | 99.809586 |
| 4 | Male | 177.499761 | 93.598695 |
| ... | ... | ... | ... |
| 195 | Male | 183.628648 | 94.562609 |
| 196 | Male | 170.823445 | 82.000329 |
| 197 | Male | 192.898929 | 105.199296 |
| 198 | Male | 168.443229 | 77.380072 |
| 199 | Male | 180.736376 | 85.398496 |
200 rows × 3 columns
height_list = male_df['Height'].tolist()
weight_list = male_df['Weight'].tolist()
height_list
[187.57142322448811, 174.70603627656138, 188.23966769513365, 182.19668514447775, 177.49976148723286, 170.8226598470285, 174.71410637922045, 173.60522939943328, 170.22813214372283, 161.17949470505067, 180.83627099875525, 181.96764500282032, 164.50647599885, 175.9789980456698, 175.87907986372045, 171.81987403906984, 183.942524230009, 162.49478733149493, 176.88575214537084, 172.55745232918412, 172.50422749230435, 176.37617726482034, 168.01879518063035, 191.02317298012184, 172.4487187000549, 173.08584330743003, 175.48766382947102, 184.91414253191076, 171.24995539556042, 173.98089500648203, 174.29000079752274, 188.04587135903623, 181.6816687303501, 175.71760919290577, 176.7260940081211, 178.81835937672955, 175.45347476554045, 170.6715491381159, 167.15058176352375, 163.3326372445611, 172.65481401707962, 183.36114193640202, 165.7866776310601, 167.86905053573216, 171.47621664991306, 178.06615708809903, 173.35966455980267, 183.31868115412516, 175.7168383519486, 185.0907147739769, 164.54776076408726, 178.26621671420858, 173.96828364246608, 171.02031043778115, 170.15699796960894, 168.90885674663477, 173.61676388185296, 179.7671260388094, 180.88844612045608, 177.8339479866312, 181.3640828407276, 176.66209294087304, 187.49907777673553, 170.172082190715, 181.40290336845038, 165.80942253087093, 173.41644634454013, 184.8240396637097, 172.97242469288085, 174.99763342330976, 175.99535530679523, 179.12899028957952, 176.94586207535545, 171.6919827721059, 184.17062576546343, 180.98258343367382, 182.39533469899166, 183.5027203172998, 168.94207425877414, 167.71374637223823, 171.6459710145469, 173.34628424812547, 187.51836258668223, 177.29825947564564, 167.53732757364295, 180.54070464239254, 169.15433920450758, 175.13689913577448, 173.34093295245825, 183.68473581202682, 182.41116794180826, 165.70210577785005, 179.42694643045218, 164.4190910170277, 170.44302000331504, 165.39841161563606, 182.12113441846498, 169.7555106430556, 168.83703916425623, 176.3052766003339, 177.93253078529222, 169.52875581598627, 168.33679595315257, 173.5542579020041, 177.99275601589463, 174.5817898977854, 171.59237020516503, 168.284226638376, 175.72279480866203, 171.72712655176483, 176.00355457559223, 173.20433737956765, 181.88199792918212, 175.760018076242, 170.84538088681722, 187.1603912359995, 176.62452611852882, 173.51136199732946, 172.05696720475584, 182.13263700679425, 161.62595209361822, 174.55184407744372, 170.0513277914747, 159.27575392770365, 183.53635866594223, 180.1079787791757, 171.72900274827228, 179.83595974746328, 176.03411337448668, 168.25696178653445, 171.43017016227572, 167.14787066005033, 181.4670123803187, 173.89959078867784, 162.5293284966492, 180.34482012204785, 173.07489944908338, 173.72431663022735, 172.96439815234191, 173.0770995040391, 174.90462938420796, 168.01808545726598, 168.16332141451562, 171.27758534530702, 176.4566093193674, 179.11028362363936, 181.19947514148774, 180.34493970907866, 168.15396317393342, 184.2599998952873, 171.3975534654244, 165.99004281175817, 179.94391854905604, 177.64372715865045, 164.43913132680018, 176.04334331619975, 174.84969566407588, 157.32080109991458, 174.22707225919999, 165.6551819683161, 163.40746619174064, 174.6222230224191, 184.12372471703657, 170.77418453236794, 170.84955630188296, 165.40100087926066, 168.3075628969345, 171.96226633545336, 166.44796775540533, 177.46718331537843, 173.94396342642182, 183.42346689812754, 173.17601304294362, 182.83226461889106, 167.80544730113522, 169.30665699133897, 171.7025374560371, 173.34469799833073, 164.6138984594212, 171.43023440539514, 173.17905599543812, 176.66559765337107, 168.66171415763515, 169.14407907130675, 182.71715192043442, 173.44908598881977, 177.68684603571606, 182.50186802932285, 167.01896654731118, 170.2609042806163, 194.8429805389624, 184.33064828309134, 177.11116443931243, 184.91002897164842, 184.24995934273218, 183.62864819641894, 170.82344503229407, 192.89892935504406, 168.4432290736113, 180.73637586471077]
weight_list
[109.72107461075917, 73.62279190675635, 96.49762886770051, 99.80958560542733, 93.59869511400913, 69.04227247266442, 83.42828690099842, 76.19041409846824, 79.80025182186384, 70.94169987966318, 84.64257043859637, 96.95136363527362, 75.80774116040074, 85.93134240266038, 84.56511612173202, 78.10267768808554, 88.91703477400334, 78.41862298582492, 84.36090409757549, 82.74733562347498, 78.97765692678098, 89.68946415369008, 67.66399134690883, 103.7645982346668, 73.48498116762939, 87.24576019177793, 83.65838787127223, 93.8156886219873, 79.47569831806146, 70.00864338247818, 85.05167337826077, 96.5744537026371, 88.46513743485168, 93.06972508941506, 92.60728954731077, 87.49956258404077, 89.57916000758136, 83.37525498065008, 74.32193757551913, 73.98455206397283, 78.07939358912762, 88.01751509532019, 76.48372312486332, 73.11611111993703, 74.68864540903805, 85.69371523224486, 84.84923918948712, 94.83294842313065, 87.09623747350628, 95.86334401413797, 75.12017001956083, 91.20468787122894, 78.66378518840078, 82.28513918350808, 76.99172897804684, 74.07595598956712, 86.05110374099618, 87.14636484363656, 89.90131106492622, 95.0396216389266, 90.15592541329404, 89.84736058288759, 88.58233162188917, 74.8170237220387, 81.58489373953027, 70.42040626611075, 83.69658798739795, 100.09787583519599, 83.14922177179501, 89.10881633176783, 83.7311594605125, 94.13531207645366, 80.37698923897868, 73.98454528811837, 98.07493943404084, 92.83017047638785, 91.13603808120453, 100.19935733758183, 89.10815773814433, 76.49418470431661, 82.29649428888571, 90.11005845518574, 107.9172160828461, 78.4901933756455, 72.88499856956516, 85.54884505018158, 94.50401713667996, 87.74069381999635, 78.97482864646162, 89.52488616116239, 91.4535960710885, 82.10564994047328, 82.60937004535216, 80.53499118796117, 74.83121104152005, 75.16801766786654, 87.58603927466554, 81.95683251192315, 78.36884245682175, 80.49931742100733, 91.08984699832631, 76.0854822416895, 70.03267736124303, 80.7325148712439, 83.43326206211684, 81.28592181382977, 71.92522920428186, 75.52473153129206, 78.01048392932307, 79.09026272821687, 85.15985910819172, 91.39340238269175, 92.95332911215124, 81.05040050815454, 74.3956929936653, 101.27357045131782, 92.06033915212032, 85.90666889334867, 80.81680456905312, 101.28687469084677, 75.7629469111835, 81.85921718454048, 74.5883417836223, 68.0250746475961, 87.50262500252731, 91.57304143567029, 81.58676419352868, 89.11626086335811, 85.07199719020707, 77.87870721462966, 86.95347806218668, 75.22858839479031, 84.2152028577875, 80.86726452093582, 68.48294567689975, 95.22233934733885, 78.99665122525411, 81.59173153702342, 76.08882335923862, 83.02766131402728, 76.90968018824113, 78.66728343510957, 82.05012178353034, 83.23825176974468, 93.43000843969563, 83.35206369176677, 91.60501823088019, 98.03042597634044, 72.29823470219306, 95.40699735969999, 81.01992077454297, 66.56005258294383, 88.87553499292459, 77.19046884781403, 76.95378448893794, 89.98306180729176, 81.60476769720982, 66.79768457517916, 80.27870591877974, 84.1267736868605, 80.50983446320866, 77.03244797388919, 90.07514314775118, 92.23816233198148, 87.6420620801761, 80.80065196349148, 73.835577024146, 72.25114585157111, 78.06973866887638, 91.55174250440393, 77.27051689588949, 93.16322833356082, 77.01443611857734, 90.91683471915486, 81.00658691755643, 76.20988194225913, 80.63349440317036, 87.66925553792113, 73.65412067823293, 83.18150352538525, 82.12566336343181, 95.64983905511497, 75.10060789186987, 88.97144023757716, 91.64671683468656, 77.29830983656412, 91.60211062665059, 87.718860650477, 81.5811685105543, 70.9826404969754, 106.61027265651639, 104.58049948085723, 92.48551337796762, 98.9812925243504, 101.0001329546353, 94.56260850473717, 82.00032886452549, 105.19929633790117, 77.38007240504437, 85.39849606118273]
xs = np.array(height_list, dtype = np.float128)
ys = np.array(weight_list, dtype = np.float128)
xs
array([187.57142322, 174.70603628, 188.2396677 , 182.19668514,
177.49976149, 170.82265985, 174.71410638, 173.6052294 ,
170.22813214, 161.17949471, 180.836271 , 181.967645 ,
164.506476 , 175.97899805, 175.87907986, 171.81987404,
183.94252423, 162.49478733, 176.88575215, 172.55745233,
172.50422749, 176.37617726, 168.01879518, 191.02317298,
172.4487187 , 173.08584331, 175.48766383, 184.91414253,
171.2499554 , 173.98089501, 174.2900008 , 188.04587136,
181.68166873, 175.71760919, 176.72609401, 178.81835938,
175.45347477, 170.67154914, 167.15058176, 163.33263724,
172.65481402, 183.36114194, 165.78667763, 167.86905054,
171.47621665, 178.06615709, 173.35966456, 183.31868115,
175.71683835, 185.09071477, 164.54776076, 178.26621671,
173.96828364, 171.02031044, 170.15699797, 168.90885675,
173.61676388, 179.76712604, 180.88844612, 177.83394799,
181.36408284, 176.66209294, 187.49907778, 170.17208219,
181.40290337, 165.80942253, 173.41644634, 184.82403966,
172.97242469, 174.99763342, 175.99535531, 179.12899029,
176.94586208, 171.69198277, 184.17062577, 180.98258343,
182.3953347 , 183.50272032, 168.94207426, 167.71374637,
171.64597101, 173.34628425, 187.51836259, 177.29825948,
167.53732757, 180.54070464, 169.1543392 , 175.13689914,
173.34093295, 183.68473581, 182.41116794, 165.70210578,
179.42694643, 164.41909102, 170.44302 , 165.39841162,
182.12113442, 169.75551064, 168.83703916, 176.3052766 ,
177.93253079, 169.52875582, 168.33679595, 173.5542579 ,
177.99275602, 174.5817899 , 171.59237021, 168.28422664,
175.72279481, 171.72712655, 176.00355458, 173.20433738,
181.88199793, 175.76001808, 170.84538089, 187.16039124,
176.62452612, 173.511362 , 172.0569672 , 182.13263701,
161.62595209, 174.55184408, 170.05132779, 159.27575393,
183.53635867, 180.10797878, 171.72900275, 179.83595975,
176.03411337, 168.25696179, 171.43017016, 167.14787066,
181.46701238, 173.89959079, 162.5293285 , 180.34482012,
173.07489945, 173.72431663, 172.96439815, 173.0770995 ,
174.90462938, 168.01808546, 168.16332141, 171.27758535,
176.45660932, 179.11028362, 181.19947514, 180.34493971,
168.15396317, 184.2599999 , 171.39755347, 165.99004281,
179.94391855, 177.64372716, 164.43913133, 176.04334332,
174.84969566, 157.3208011 , 174.22707226, 165.65518197,
163.40746619, 174.62222302, 184.12372472, 170.77418453,
170.8495563 , 165.40100088, 168.3075629 , 171.96226634,
166.44796776, 177.46718332, 173.94396343, 183.4234669 ,
173.17601304, 182.83226462, 167.8054473 , 169.30665699,
171.70253746, 173.344698 , 164.61389846, 171.43023441,
173.179056 , 176.66559765, 168.66171416, 169.14407907,
182.71715192, 173.44908599, 177.68684604, 182.50186803,
167.01896655, 170.26090428, 194.84298054, 184.33064828,
177.11116444, 184.91002897, 184.24995934, 183.6286482 ,
170.82344503, 192.89892936, 168.44322907, 180.73637586],
dtype=float128)
ys
array([109.72107461, 73.62279191, 96.49762887, 99.80958561,
93.59869511, 69.04227247, 83.4282869 , 76.1904141 ,
79.80025182, 70.94169988, 84.64257044, 96.95136364,
75.80774116, 85.9313424 , 84.56511612, 78.10267769,
88.91703477, 78.41862299, 84.3609041 , 82.74733562,
78.97765693, 89.68946415, 67.66399135, 103.76459823,
73.48498117, 87.24576019, 83.65838787, 93.81568862,
79.47569832, 70.00864338, 85.05167338, 96.5744537 ,
88.46513743, 93.06972509, 92.60728955, 87.49956258,
89.57916001, 83.37525498, 74.32193758, 73.98455206,
78.07939359, 88.0175151 , 76.48372312, 73.11611112,
74.68864541, 85.69371523, 84.84923919, 94.83294842,
87.09623747, 95.86334401, 75.12017002, 91.20468787,
78.66378519, 82.28513918, 76.99172898, 74.07595599,
86.05110374, 87.14636484, 89.90131106, 95.03962164,
90.15592541, 89.84736058, 88.58233162, 74.81702372,
81.58489374, 70.42040627, 83.69658799, 100.09787584,
83.14922177, 89.10881633, 83.73115946, 94.13531208,
80.37698924, 73.98454529, 98.07493943, 92.83017048,
91.13603808, 100.19935734, 89.10815774, 76.4941847 ,
82.29649429, 90.11005846, 107.91721608, 78.49019338,
72.88499857, 85.54884505, 94.50401714, 87.74069382,
78.97482865, 89.52488616, 91.45359607, 82.10564994,
82.60937005, 80.53499119, 74.83121104, 75.16801767,
87.58603927, 81.95683251, 78.36884246, 80.49931742,
91.089847 , 76.08548224, 70.03267736, 80.73251487,
83.43326206, 81.28592181, 71.9252292 , 75.52473153,
78.01048393, 79.09026273, 85.15985911, 91.39340238,
92.95332911, 81.05040051, 74.39569299, 101.27357045,
92.06033915, 85.90666889, 80.81680457, 101.28687469,
75.76294691, 81.85921718, 74.58834178, 68.02507465,
87.502625 , 91.57304144, 81.58676419, 89.11626086,
85.07199719, 77.87870721, 86.95347806, 75.22858839,
84.21520286, 80.86726452, 68.48294568, 95.22233935,
78.99665123, 81.59173154, 76.08882336, 83.02766131,
76.90968019, 78.66728344, 82.05012178, 83.23825177,
93.43000844, 83.35206369, 91.60501823, 98.03042598,
72.2982347 , 95.40699736, 81.01992077, 66.56005258,
88.87553499, 77.19046885, 76.95378449, 89.98306181,
81.6047677 , 66.79768458, 80.27870592, 84.12677369,
80.50983446, 77.03244797, 90.07514315, 92.23816233,
87.64206208, 80.80065196, 73.83557702, 72.25114585,
78.06973867, 91.5517425 , 77.2705169 , 93.16322833,
77.01443612, 90.91683472, 81.00658692, 76.20988194,
80.6334944 , 87.66925554, 73.65412068, 83.18150353,
82.12566336, 95.64983906, 75.10060789, 88.97144024,
91.64671683, 77.29830984, 91.60211063, 87.71886065,
81.58116851, 70.9826405 , 106.61027266, 104.58049948,
92.48551338, 98.98129252, 101.00013295, 94.5626085 ,
82.00032886, 105.19929634, 77.38007241, 85.39849606],
dtype=float128)
style.use('seaborn')
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
<matplotlib.collections.PathCollection at 0x7fcb12e70670>
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
plt.show()
def best_fit_line(xs, ys):
a1 = (mean(xs) * mean(ys) - mean(xs * ys))/(mean(xs) * mean(xs) - mean(xs * xs))
a0 = mean(ys) - a1 * mean(xs)
return a1, a0
slope, y_intercept = best_fit_line(xs, ys)
regression_line = [slope * x + y_intercept for x in xs]
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
plt.plot(xs, regression_line, label = "Best Fit Line", color = 'orange', linewidth = 4)
plt.show()
average_man_heghit = 170
average_man_weight = slope * average_man_heghit + y_intercept
average_man_weight
79.23617231210749291
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
plt.plot(xs, regression_line, label = "Best Fit Line", color = 'orange', linewidth = 4)
plt.scatter(average_man_heghit, average_man_weight, label = "Average Man Prediction", color = 'red', s = 100)
plt.title("Height and Wieight Linear Regression")
plt.xlabel('Height (cm)')
plt.ylabel('Weight (kg)')
plt.legend()
plt.show()
from sklearn import linear_model
from sklearn.linear_model import LinearRegression
height_weight = linear_model.LinearRegression()
height_weight.fit(xs.reshape(-1, 1), ys)
LinearRegression()
regression_line = height_weight.predict(xs.reshape(-1, 1))
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
plt.plot(xs, regression_line, label = "Best Fit Line", color = 'orange', linewidth = 4)
# plt.scatter(average_man_heghit, average_man_weight, label = "Average Man Prediction", color = 'red', s = 100)
plt.title("Height and Wieight Linear Regression")
plt.xlabel('Height (cm)')
plt.ylabel('Weight (kg)')
plt.legend()
plt.show()
KIS_height = 170
KIS_wieght = height_weight.predict(np.array([[KIS_height]]))[0]
KIS_wieght
79.23617231210784
plt.scatter(xs, ys, label = "Data Points", alpha = 0.6, color = 'green', s = 75)
plt.plot(xs, regression_line, label = "Best Fit Line", color = 'orange', linewidth = 4)
plt.scatter(KIS_height, KIS_wieght, label = "KSI Man Prediction", color = 'red', s = 100)
plt.title("Height and Wieight Linear Regression")
plt.xlabel('Height (cm)')
plt.ylabel('Weight (kg)')
plt.legend()
plt.show()
import numpy as np
x = np.array([1, 2, 3, 4, 5])
x
array([1, 2, 3, 4, 5])
from numpy import *
x = array([0, 1, 2, 3, 4, 5])
x
array([0, 1, 2, 3, 4, 5])
y = array([0, 0.8, 0.9, 0.1, -0.8, -1])
y
array([ 0. , 0.8, 0.9, 0.1, -0.8, -1. ])
from scipy.interpolate import *
p1 = polyfit(x, y, 1)
print(p1)
[-0.30285714 0.75714286]
from matplotlib.pylab import *
plot(x, y, "o")
show()
plot(x, polyval(p1, x))
[<matplotlib.lines.Line2D at 0x7fcb08d652b0>]
plot(x, y, "o")
plot(x, polyval(p1, x), 'r-')
[<matplotlib.lines.Line2D at 0x7fcb130baf40>]
p2 = polyfit(x, y, 2)
print(p1)
[-0.30285714 0.75714286]
print(p1[0])
-0.3028571428571428
print(p1[1])
0.7571428571428572
print(p2)
[-0.16071429 0.50071429 0.22142857]
print(p2[0], p2[1], p2[2])
-0.16071428571428573 0.5007142857142854 0.22142857142857117
plot(x, y, "o")
plot(x, polyval(p1, x), 'r-')
plot(x, polyval(p2, x), 'b--')
[<matplotlib.lines.Line2D at 0x7fcb13604e50>]
p3 = polyfit(x, y, 3)
print(p1)
[-0.30285714 0.75714286]
print(p2)
[-0.16071429 0.50071429 0.22142857]
print(p3)
[ 0.08703704 -0.81349206 1.69312169 -0.03968254]
print(p3[0], p3[1], p3[2], p3[3])
0.08703703703703657 -0.8134920634920602 1.6931216931216868 -0.03968253968253445
plot(x, y, "o")
plot(x, polyval(p1, x), 'r-')
plot(x, polyval(p2, x), 'b--')
plot(x, polyval(p3, x), 'm--')
show()
plot(x, y, "o")
xp = linspace(-2, 6, 100)
plot(xp, polyval(p1, xp), 'r-')
plot(xp, polyval(p2, xp), 'b--')
plot(xp, polyval(p3, xp), 'm--')
show()
yfit = p1[0] * x + p1[1]
print(yfit)
print(y)
[ 0.75714286 0.45428571 0.15142857 -0.15142857 -0.45428571 -0.75714286] [ 0. 0.8 0.9 0.1 -0.8 -1. ]
scatter(x, y)
plot(x, yfit, 'r')
[<matplotlib.lines.Line2D at 0x7fcb137f4250>]
SS_res = sum(pow(y - yfit, 2))
SS_tot = sum(pow(y - mean(y), 2))
R_Squared = 1 - SS_res/SS_tot
print(R_Squared)
0.5177880184331798
from scipy.stats import *
slope, intercept, r_value, p_value, std_err = linregress(x, y)
print(pow(r_value, 2))
0.5177880184331796
print(slope)
-0.3028571428571429
print(p1[0])
-0.3028571428571428
print(intercept)
0.7571428571428572
print(p1[1])
0.7571428571428572
print(std_err)
0.14613384974336482
print(p_value)
0.10693129655435198
a = np.array([1,2,3,4,5])
print(a)
[1 2 3 4 5]
s = np.sum(a)
print(s)
15
avg = np.mean(a)
print(avg)
3.0
np.std(a)
1.4142135623730951
b = a - avg
print(b)
[-2. -1. 0. 1. 2.]
c = b ** 2
print(c)
[4. 1. 0. 1. 4.]
d = sum(c)
print(d)
10.0
e = d / 5
print(e)
2.0
f = np.sqrt(e)
print(f)
1.4142135623730951
ages = [5,31,43,48,50,41,7,11,15,39,80,82,32,2,8,6,25,36,27,61,31]
x = np.percentile(ages, 75)
print(x)
43.0
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import style
import scipy as sp
from scipy.interpolate import interp1d
print(plt.style.available)
['Solarize_Light2', '_classic_test_patch', 'bmh', 'bright', 'cjk-jp-font', 'cjk-kr-font', 'cjk-sc-font', 'cjk-tc-font', 'classic', 'dark_background', 'fast', 'fivethirtyeight', 'ggplot', 'grayscale', 'grid', 'high-contrast', 'high-vis', 'ieee', 'latex-sans', 'light', 'muted', 'nature', 'no-latex', 'notebook', 'pgf', 'retro', 'scatter', 'science', 'seaborn', 'seaborn-bright', 'seaborn-colorblind', 'seaborn-dark', 'seaborn-dark-palette', 'seaborn-darkgrid', 'seaborn-deep', 'seaborn-muted', 'seaborn-notebook', 'seaborn-paper', 'seaborn-pastel', 'seaborn-poster', 'seaborn-talk', 'seaborn-ticks', 'seaborn-white', 'seaborn-whitegrid', 'std-colors', 'tableau-colorblind10', 'vibrant']
plt.style.use(['seaborn', 'notebook', 'grid'])
plt.figure(figsize=(10,4))
x = np.linspace(0, 10, 10)
y = x**2 * np.sin(x)
plt.scatter(x,y, color = 'green', alpha = 0.6, s = 85, label='Experimental Data')
plt.legend(loc='lower left', fontsize=12)
show()
f = interp1d(x, y, kind='cubic')
x_dense = np.linspace(0, 10, 100)
y_dense = f(x_dense)
plt.figure(figsize=(10,4))
plt.plot(x_dense, y_dense, '-', color = 'orange', alpha = 0.8, lw=2, zorder = 0, label='Fitted Curve')
plt.scatter(x,y, color = 'green', alpha = 0.6, s = 85, zorder = 1, label='Experimental Data')
plt.xlabel('$\sigma^2$(dimenssionless)', fontsize=16)
plt.ylabel('$\int_0^\infty f(x)dx$', fontsize=16)
plt.title('Title', fontsize=20)
plt.legend(loc='lower left', fontsize=12)
show()
x_data = np.linspace(0, 10, 10)
y_data = 3*x_data**2 + 2
plt.figure(figsize=(10,4))
plt.scatter(x_data, y_data, color = 'green', alpha = 0.6, s = 85, zorder = 1, label='Experimental Data')
<matplotlib.collections.PathCollection at 0x7fcad1ac7c10>
from scipy.optimize import curve_fit
def func(x, a, b):
return a*x**2 + b
popt, pcov = curve_fit(func, x_data, y_data, p0=(1,1))
popt
array([3., 2.])
t_data = np.array([ 0. , 0.34482759, 0.68965517, 1.03448276, 1.37931034,
1.72413793, 2.06896552, 2.4137931 , 2.75862069, 3.10344828,
3.44827586, 3.79310345, 4.13793103, 4.48275862, 4.82758621,
5.17241379, 5.51724138, 5.86206897, 6.20689655, 6.55172414,
6.89655172, 7.24137931, 7.5862069 , 7.93103448, 8.27586207,
8.62068966, 8.96551724, 9.31034483, 9.65517241, 10. ])
y_data = np.array([ 4.3303953 , 1.61137995, -2.15418696, -3.90137249, -1.67259042,
2.16884383, 3.86635998, 1.85194506, -1.8489224 , -3.96560495,
-2.13385255, 1.59425817, 4.06145238, 1.89300594, -1.76870297,
-4.26791226, -2.46874133, 1.37019912, 4.24945607, 2.27038039,
-1.50299303, -3.46774049, -2.50845488, 1.20022052, 3.81633703,
2.91511556, -1.24569189, -3.72716214, -2.54549857, 0.87262548])
plt.plot(t_data,y_data,'o--', color = 'red', alpha = 0.6)
[<matplotlib.lines.Line2D at 0x7fcad1bedf40>]
from scipy.optimize import curve_fit
def func(x, A, w, phi):
return A*np.cos(w*x+phi)
popt, pcov = curve_fit(func, t_data, y_data, p0=(4, np.pi, 0))
popt
array([3.94836218, 2.99899521, 0.10411349])
A, w, phi = popt
t = np.linspace(0, 10, 100)
y = func(t, A, w, phi)
plt.scatter(t_data,y_data)
plt.plot(t,y)
[<matplotlib.lines.Line2D at 0x7fcad6850f40>]
popt
array([3.94836218, 2.99899521, 0.10411349])
np.sqrt(np.diag(pcov))
array([0.05117448, 0.00430857, 0.02575703])
help(curve_fit)
Help on function curve_fit in module scipy.optimize.minpack:
curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs)
Use non-linear least squares to fit a function, f, to data.
Assumes ``ydata = f(xdata, *params) + eps``.
Parameters
----------
f : callable
The model function, f(x, ...). It must take the independent
variable as the first argument and the parameters to fit as
separate remaining arguments.
xdata : array_like or object
The independent variable where the data is measured.
Should usually be an M-length sequence or an (k,M)-shaped array for
functions with k predictors, but can actually be any object.
ydata : array_like
The dependent data, a length M array - nominally ``f(xdata, ...)``.
p0 : array_like, optional
Initial guess for the parameters (length N). If None, then the
initial values will all be 1 (if the number of parameters for the
function can be determined using introspection, otherwise a
ValueError is raised).
sigma : None or M-length sequence or MxM array, optional
Determines the uncertainty in `ydata`. If we define residuals as
``r = ydata - f(xdata, *popt)``, then the interpretation of `sigma`
depends on its number of dimensions:
- A 1-D `sigma` should contain values of standard deviations of
errors in `ydata`. In this case, the optimized function is
``chisq = sum((r / sigma) ** 2)``.
- A 2-D `sigma` should contain the covariance matrix of
errors in `ydata`. In this case, the optimized function is
``chisq = r.T @ inv(sigma) @ r``.
.. versionadded:: 0.19
None (default) is equivalent of 1-D `sigma` filled with ones.
absolute_sigma : bool, optional
If True, `sigma` is used in an absolute sense and the estimated parameter
covariance `pcov` reflects these absolute values.
If False (default), only the relative magnitudes of the `sigma` values matter.
The returned parameter covariance matrix `pcov` is based on scaling
`sigma` by a constant factor. This constant is set by demanding that the
reduced `chisq` for the optimal parameters `popt` when using the
*scaled* `sigma` equals unity. In other words, `sigma` is scaled to
match the sample variance of the residuals after the fit. Default is False.
Mathematically,
``pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)``
check_finite : bool, optional
If True, check that the input arrays do not contain nans of infs,
and raise a ValueError if they do. Setting this parameter to
False may silently produce nonsensical results if the input arrays
do contain nans. Default is True.
bounds : 2-tuple of array_like, optional
Lower and upper bounds on parameters. Defaults to no bounds.
Each element of the tuple must be either an array with the length equal
to the number of parameters, or a scalar (in which case the bound is
taken to be the same for all parameters). Use ``np.inf`` with an
appropriate sign to disable bounds on all or some parameters.
.. versionadded:: 0.17
method : {'lm', 'trf', 'dogbox'}, optional
Method to use for optimization. See `least_squares` for more details.
Default is 'lm' for unconstrained problems and 'trf' if `bounds` are
provided. The method 'lm' won't work when the number of observations
is less than the number of variables, use 'trf' or 'dogbox' in this
case.
.. versionadded:: 0.17
jac : callable, string or None, optional
Function with signature ``jac(x, ...)`` which computes the Jacobian
matrix of the model function with respect to parameters as a dense
array_like structure. It will be scaled according to provided `sigma`.
If None (default), the Jacobian will be estimated numerically.
String keywords for 'trf' and 'dogbox' methods can be used to select
a finite difference scheme, see `least_squares`.
.. versionadded:: 0.18
kwargs
Keyword arguments passed to `leastsq` for ``method='lm'`` or
`least_squares` otherwise.
Returns
-------
popt : array
Optimal values for the parameters so that the sum of the squared
residuals of ``f(xdata, *popt) - ydata`` is minimized.
pcov : 2-D array
The estimated covariance of popt. The diagonals provide the variance
of the parameter estimate. To compute one standard deviation errors
on the parameters use ``perr = np.sqrt(np.diag(pcov))``.
How the `sigma` parameter affects the estimated covariance
depends on `absolute_sigma` argument, as described above.
If the Jacobian matrix at the solution doesn't have a full rank, then
'lm' method returns a matrix filled with ``np.inf``, on the other hand
'trf' and 'dogbox' methods use Moore-Penrose pseudoinverse to compute
the covariance matrix.
Raises
------
ValueError
if either `ydata` or `xdata` contain NaNs, or if incompatible options
are used.
RuntimeError
if the least-squares minimization fails.
OptimizeWarning
if covariance of the parameters can not be estimated.
See Also
--------
least_squares : Minimize the sum of squares of nonlinear functions.
scipy.stats.linregress : Calculate a linear least squares regression for
two sets of measurements.
Notes
-----
With ``method='lm'``, the algorithm uses the Levenberg-Marquardt algorithm
through `leastsq`. Note that this algorithm can only deal with
unconstrained problems.
Box constraints can be handled by methods 'trf' and 'dogbox'. Refer to
the docstring of `least_squares` for more information.
Examples
--------
>>> import matplotlib.pyplot as plt
>>> from scipy.optimize import curve_fit
>>> def func(x, a, b, c):
... return a * np.exp(-b * x) + c
Define the data to be fit with some noise:
>>> xdata = np.linspace(0, 4, 50)
>>> y = func(xdata, 2.5, 1.3, 0.5)
>>> rng = np.random.default_rng()
>>> y_noise = 0.2 * rng.normal(size=xdata.size)
>>> ydata = y + y_noise
>>> plt.plot(xdata, ydata, 'b-', label='data')
Fit for the parameters a, b, c of the function `func`:
>>> popt, pcov = curve_fit(func, xdata, ydata)
>>> popt
array([2.56274217, 1.37268521, 0.47427475])
>>> plt.plot(xdata, func(xdata, *popt), 'r-',
... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
Constrain the optimization to the region of ``0 <= a <= 3``,
``0 <= b <= 1`` and ``0 <= c <= 0.5``:
>>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5]))
>>> popt
array([2.43736712, 1. , 0.34463856])
>>> plt.plot(xdata, func(xdata, *popt), 'g--',
... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
>>> plt.xlabel('x')
>>> plt.ylabel('y')
>>> plt.legend()
>>> plt.show()
help(diag)
Help on function diag in module numpy:
diag(v, k=0)
Extract a diagonal or construct a diagonal array.
See the more detailed documentation for ``numpy.diagonal`` if you use this
function to extract a diagonal and wish to write to the resulting array;
whether it returns a copy or a view depends on what version of numpy you
are using.
Parameters
----------
v : array_like
If `v` is a 2-D array, return a copy of its `k`-th diagonal.
If `v` is a 1-D array, return a 2-D array with `v` on the `k`-th
diagonal.
k : int, optional
Diagonal in question. The default is 0. Use `k>0` for diagonals
above the main diagonal, and `k<0` for diagonals below the main
diagonal.
Returns
-------
out : ndarray
The extracted diagonal or constructed diagonal array.
See Also
--------
diagonal : Return specified diagonals.
diagflat : Create a 2-D array with the flattened input as a diagonal.
trace : Sum along diagonals.
triu : Upper triangle of an array.
tril : Lower triangle of an array.
Examples
--------
>>> x = np.arange(9).reshape((3,3))
>>> x
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
>>> np.diag(x)
array([0, 4, 8])
>>> np.diag(x, k=1)
array([1, 5])
>>> np.diag(x, k=-1)
array([3, 7])
>>> np.diag(np.diag(x))
array([[0, 0, 0],
[0, 4, 0],
[0, 0, 8]])
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit
Given some $x_{\text{data}}=[...]$ and $y_{\text{data}}=[...]$ and a model function $f$ that depends on unknown parameters $\beta$ the goal is to find the optimal set of parameters $\beta$ such that the function $y=f(x,\beta)$ best resembles the data.
One way of doing this is to minimize $\sum_i (f(x_i,\beta) - y_i)^2$ by adjusting the parameters in $\beta$. Here $x_i$ and $y_i$ are the ith data point. This is called the method of least squares
If the $y_i$ data points have corresponding errors on them, then minimizing $\sum_i (f(x_i,\beta) - y_i)^2/\sigma_i^2$ gives the so called maximum-likelihood estimate for $\beta$. Given the data $x_i$ and $y_i$, the the value of $\beta$ obtained through this method is the most likely configuration.
x_data = np.array([0. , 0.15789474, 0.31578947, 0.47368421, 0.63157895,
0.78947368, 0.94736842, 1.10526316, 1.26315789, 1.42105263,
1.57894737, 1.73684211, 1.89473684, 2.05263158, 2.21052632,
2.36842105, 2.52631579, 2.68421053, 2.84210526, 3. ])
y_data = np.array([ 2.95258285, 2.49719803, -2.1984975 , -4.88744346,
-7.41326345, -8.44574157, -10.01878504, -13.83743553,
-12.91548145, -15.41149046, -14.93516299, -13.42514157,
-14.12110495, -17.6412464 , -16.1275509 , -16.11533771,
-15.66076021, -13.48938865, -11.33918701, -11.70467566])
plt.scatter(x_data,y_data)
plt.show()
Here is a model where $\beta = (a,b,c)$ and
$$f(x, \beta) = f(x, a, b, c) = a(x-b)^2 + c$$We want to find the optimal values of $a$, $b$, and $c$ that fit the data above
pcov: the covariance matrix, which gives an estimate of the "error" of the parameters (based on how sensitive the "goodness-of-fit" is to changing them) and the relation the parameters have to eachother (are two seperate parameters really needed?)def model_f(x, a, b, c):
return a*(x-b)**2 + c
popt, pcov = curve_fit(model_f, x_data, y_data, p0=[3,2,-16])
popt
array([ 4.92207991, 2.00492121, -16.04690266])
pcov
array([[ 0.12827822, -0.01315915, -0.07363194],
[-0.01315915, 0.00222114, 0.00322289],
[-0.07363194, 0.00322289, 0.13377697]])
a_opt, b_opt, c_opt = popt
x_model = np.linspace(min(x_data), max(x_data), 100)
y_model = model_f(x_model, a_opt, b_opt, c_opt)
plt.scatter(x_data,y_data)
plt.plot(x_model,y_model, color='r')
plt.show()
pcov is difficult to interpret, but we can do so as follows.
Lets look at an example below
def fit_f(x, a, b, c, d):
return a*(x-b)**2 + c + d*0.0001*np.cos(x)
popt, pcov = curve_fit(fit_f, x_data, y_data, p0=[1,2,-16,1])
pcov
array([[ 1.45024166e-01, -4.03732456e-02, 4.08922209e-02,
3.05314767e+03],
[-4.03732456e-02, 7.20537190e-02, -3.07488450e-01,
-8.46751020e+03],
[ 4.08922209e-02, -3.07488450e-01, 1.52776292e+00,
3.78842445e+04],
[ 3.05314767e+03, -8.46751020e+03, 3.78842445e+04,
1.03326900e+09]])
style.use('default')
plt.imshow(np.log(np.abs(pcov)))
plt.colorbar()
plt.show()
df = pd.read_csv('data/curvefite1.csv')
x_stuff = df.Channel
y_stuff = df.Cs137
plt.scatter(x_stuff, y_stuff)
<matplotlib.collections.PathCollection at 0x7fcad1a69c40>
def gauss_f(x,A,mu,sig):
return A*np.exp(-(x-mu)**2/sig**2)
x_data = x_stuff[(x_stuff>725) & (x_stuff<850)]
y_data = y_stuff[(x_stuff>725) & (x_stuff<850)]
plt.scatter(x_data, y_data)
plt.show()
popt, pcov = curve_fit(gauss_f, x_data, y_data, p0=[2000, 790, 60])
popt
array([1939.55510368, 786.33851861, 57.78717219])
A_opt, mu_opt, sig_opt = popt
x_model = np.linspace(min(x_data), max(x_data), 1000)
y_model = gauss_f(x_model, A_opt, mu_opt, sig_opt)
plt.scatter(x_data, y_data)
plt.plot(x_model, y_model, color='r')
plt.show()
popt
array([1939.55510368, 786.33851861, 57.78717219])
np.sqrt(np.diag(pcov))
array([5.70257762, 0.14401935, 0.26427165])
style.use('default')
plt.imshow(np.log(np.abs(pcov)))
plt.colorbar()
plt.show()
pcov
array([[ 3.25193916e+01, 1.08041881e-02, -9.64139420e-01],
[ 1.08041881e-02, 2.07415740e-02, -1.14109432e-03],
[-9.64139420e-01, -1.14109432e-03, 6.98395035e-02]])
np.sqrt(np.diag(pcov))
array([5.70257762, 0.14401935, 0.26427165])
yerr_data = np.sqrt(y_data)
plt.errorbar(x_data, y_data, yerr_data, ls='', color='k')
plt.scatter(x_data, y_data, s=7, zorder=1000)
plt.show()
popt, pcov = curve_fit(gauss_f, x_data, y_data, sigma=yerr_data, p0=[2000, 790, 60])
popt
array([1936.83040984, 786.29515455, 57.88161227])
popt, pcov = curve_fit(gauss_f, x_data, y_data, p0=[2000, 790, 60])
popt
array([1939.55510368, 786.33851861, 57.78717219])
r_data, V_data = np.loadtxt('data/curvefite2.txt')
plt.scatter(r_data, V_data)
plt.xlabel('$r$ [m]', fontsize=20)
plt.ylabel('$V(r)$ [J]', fontsize=20)
plt.show()
First its probably a good idea to standardize our units. We can convert back after:
Multiply x-axis by 1e9. This converts all quantities with dimenions of meters in this problem to nanometers
Divide y-axis by Boltzmann constant $k_B=1.3806452 × 10-23J/K$ This converts all quantites with dimensions of joules to Kelvin
r_data*=1e9
V_data/=1.38e-23
plt.scatter(r_data, V_data)
plt.xlabel('$r$ [nm]', fontsize=20)
plt.ylabel('$V(r)$ [K]', fontsize=20)
plt.show()
def model_f(r, sigma, epsilon):
return 4*epsilon*((sigma/r)**12-(sigma/r)**6)
In order to fit this data to the curve, we need to identify some key points for out initial parameter guess. From wikipedia:
eps_init_guess = 120
sig_init_guess = 2**(-1/6) * 0.4
popt, pcov = curve_fit(model_f, r_data, V_data, p0=[sig_init_guess,
eps_init_guess])
sig_opt, eps_opt = popt
r_model = np.linspace(min(r_data), max(r_data), 1000)
V_model = model_f(r_model, sig_opt, eps_opt)
plt.scatter(r_data, V_data)
plt.plot(r_model, V_model, color='r')
plt.xlabel('$r$ [nm]', fontsize=20)
plt.ylabel('$V(r)$ [K]', fontsize=20)
plt.show()
sig_opt
0.3414387956041442
eps_opt
124.98902610819769
pcov
array([[ 6.41307906e-07, -2.96473497e-03],
[-2.96473497e-03, 3.49722167e+01]])
style.use('default')
plt.imshow(np.log(np.abs(pcov)))
plt.colorbar()
plt.show()