%matplotlib inline
import matplotlib.pyplot as plt


"""
Explanation: 计算传播与统计初步


分析天涯论坛的帖子


王成军 
wangchengjun@nju.edu.cn
计算传播网 http://computational-communication.com
使用pandas清洗泰坦尼克数据
分析天涯回帖数据
End of explanation
"""


import pandas as pd

"""
Explanation: 
End of explanation
"""


# Import the Pandas library
import pandas as pd
# Load the train and test datasets to create two DataFrames
train_url = "http://s3.amazonaws.com/assets.datacamp.com/course/Kaggle/train.csv"
train = pd.read_csv(train_url)

test_url = "http://s3.amazonaws.com/assets.datacamp.com/course/Kaggle/test.csv"
test = pd.read_csv(test_url)
#Print the `head` of the train and test dataframes
print(train.head())
print(test.head())

train.to_csv('/Users/chengjun/github/cjc2016/data/tatanic_train.csv')
test.to_csv('/Users/chengjun/github/cjc2016/data/tatanic_test.csv')

"""
Explanation: 使用pandas清洗泰坦尼克数据
End of explanation
"""


import pandas as pd

train = pd.read_csv('/Users/chengjun/github/cjc2016/data/tatanic_train.csv', sep = ",", header=0)

train.head()

train.describe()

train.shape

# Passengers that survived vs passengers that passed away
train["Survived"].value_counts()

# As proportions
train["Survived"].value_counts(normalize = True)

train['Sex'].value_counts()

# Males that survived vs males that passed away
train["Survived"][train["Sex"] == 'male'].value_counts()

# Females that survived vs Females that passed away
train["Survived"][train["Sex"] == 'female'].value_counts()

# Normalized male survival
train["Survived"][train["Sex"] == 'male'].value_counts(normalize = True) 

# Normalized female survival
train["Survived"][train["Sex"] == 'female'].value_counts(normalize = True)

# Create the column Child, and indicate whether child or not a child. Print the new column.
train["Child"] = float('NaN')
train.Child[train.Age < 18] = 1
train.Child[train.Age >= 18] = 0
print train.Child[:3]

# Normalized Survival Rates for under 18
train.Survived[train.Child == 1].value_counts(normalize = True)

# Normalized Survival Rates for over 18
train.Survived[train.Child == 0].value_counts(normalize = True)

# Create a copy of test: test_one
test_one = test
# Initialize a Survived column to 0
test_one['Survived'] = 0
# Set Survived to 1 if Sex equals "female" and print the `Survived` column from `test_one`
test_one.Survived[test_one.Sex =='female'] = 1

print test_one.Survived[:3]

#Convert the male and female groups to integer form
train["Sex"][train["Sex"] == "male"] = 0
train["Sex"][train["Sex"] == "female"] = 1

#Impute the Embarked variable
train["Embarked"] = train["Embarked"].fillna('S')

#Convert the Embarked classes to integer form
train["Embarked"][train["Embarked"] == "S"] = 0
train["Embarked"][train["Embarked"] == "C"] = 1
train["Embarked"][train["Embarked"] == "Q"] = 2

"""
Explanation: You can easily explore a DataFrame
 - .describe() summarizes the columns/features of the DataFrame, including the count of observations, mean, max and so on. 
 - Another useful trick is to look at the dimensions of the DataFrame. This is done by requesting the .shape attribute of your DataFrame object. (ex. your_data.shape)
从本机读取数据
End of explanation
"""


df = pd.read_csv('/Users/chengjun/github/cjc2016/data/tianya_bbs_threads_list.txt', sep = "\t", header=None)
df[:2]

df=df.rename(columns = {0:'title', 1:'link', 2:'author',3:'author_page', 4:'click', 5:'reply', 6:'time'})
df[:5]

da = pd.read_csv('/Users/chengjun/github/cjc2016/data/tianya_bbs_threads_author_info.txt', sep = "\t", header=None)
da[:2]

da=da.rename(columns = {0:'author_page', 1:'followed_num', 2:'fans_num',3:'post_num', 4:'comment_num'})
da[:5]

data = pd.concat([df,da], axis=1)
len(data)

data[:3]

"""
Explanation: 分析天涯回帖数据
End of explanation
"""


type(data.time[0])

# extract date from datetime
date = map(lambda x: x[:10], data.time)
data['date'] = pd.to_datetime(date)

# convert str to datetime format
data.time = pd.to_datetime(data.time)
data['month'] = data.time.dt.month
data['year'] = data.time.dt.year
data['day'] = data.time.dt.day
type(data.time[0])

data[:3]

"""
Explanation: Time
End of explanation
"""


import statsmodels.api as sm

'   '.join(dir(sm))  

'   '.join(dir(sm.stats))

"""
Explanation: Statsmodels
http://statsmodels.sourceforge.net/

Statsmodels is a Python module that allows users to explore data, estimate statistical models, and perform statistical tests. 
An extensive list of descriptive statistics, statistical tests, plotting functions, and result statistics are available for different types of data and each estimator. 
Researchers across fields may find that statsmodels fully meets their needs for statistical computing and data analysis in Python. 
Features include:

Linear regression models
Generalized linear models
Discrete choice models
Robust linear models
Many models and functions for time series analysis
Nonparametric estimators
A collection of datasets for examples
A wide range of statistical tests
Input-output tools for producing tables in a number of formats and for reading Stata files into NumPy and Pandas.
Plotting functions
Extensive unit tests to ensure correctness of results
Many more models and extensions in development
End of explanation
"""


data.describe()

import numpy as np

np.mean(data.click), np.std(data.click), np.sum(data.click)

# 不加权的变量描述
d1 = sm.stats.DescrStatsW(data.click, weights=[1 for i in data.click])
d1.mean, d1.var, d1.std, d1.sum

# 加权的变量描述
d1 = sm.stats.DescrStatsW(data.click, weights=data.reply)
d1.mean, d1.var, d1.std, d1.sum

np.median(data.click)

plt.hist(data.click)
plt.show()

plt.hist(data.reply, color = 'green')
plt.show()

plt.hist(np.log(data.click+1), color='green')
plt.hist(np.log(data.reply+1), color='red')
plt.show()

# Plot the height and weight to see
plt.boxplot([np.log(data.click+1)])
plt.show()

# Plot the height and weight to see
plt.boxplot([data.click, data.reply])
plt.show()

def transformData(dat):
    results = []
    for i in dat:
        if i != 'na':
            results.append( int(i))
        else:
            results.append(0)
    return results

data.fans_num = transformData(data.fans_num)
data.followed_num = transformData(data.followed_num )
data.post_num = transformData(data.post_num )
data.comment_num = transformData(data.comment_num )

# Plot the height and weight to see
plt.boxplot([np.log(data.click+1), np.log(data.reply+1), 
             np.log(data.fans_num+1), np.log(data.followed_num + 1)], 
            labels = ['$Click$', '$Reply$', '$Fans$', '$Followed$'])
plt.show()

"""
Explanation: Describe
End of explanation
"""


fig = plt.figure(figsize=(12,4))
data.boxplot(return_type='dict')
plt.yscale('log')
plt.show()

'   '.join(dir(data))

from pandas.tools import plotting

#fig = plt.figure(figsize=(10, 10))
plotting.scatter_matrix(data[['click', 'reply', 'post_num','comment_num']]) 
plt.show()

'  '.join(dir(plotting))

"""
Explanation: Pandas自身已经包含了boxplot的功能
End of explanation
"""


import seaborn # conda install seaborn
seaborn.pairplot(data, vars=['click', 'reply', 'post_num', 'comment_num'],
                  kind='reg')  

seaborn.pairplot(data, vars=['click', 'reply', 'post_num'],
                 kind='reg', hue='year')  

seaborn.lmplot(y='reply', x='click', data=data)  


"""
Explanation: 更多使用pandas.plotting绘图的操作见：
http://pandas.pydata.org/pandas-docs/version/0.15.0/visualization.html
End of explanation
"""


data.year.value_counts()

d = data.year.value_counts()
dd = pd.DataFrame(d)
dd = dd.sort_index(axis=0, ascending=True)
dd

dd.index

dd_date_str = map(lambda x: str(x) +'-01-01', dd.index)
dd_date_str

dd_date = pd.to_datetime(dd_date_str)
dd_date

plt.plot(dd_date, dd.year, 'r-o')
plt.show()

ds = dd.cumsum()
ds

d = data.year.value_counts()
dd = pd.DataFrame(d)
dd = dd.sort_index(axis=0, ascending=True)
ds = dd.cumsum()

def getDate(dat):
    dat_date_str = map(lambda x: str(x) +'-01-01', dat.index)
    dat_date = pd.to_datetime(dat_date_str)
    return dat_date

ds.date = getDate(ds)
dd.date = getDate(dd)

plt.plot(ds_date, ds.year, 'g-s', label = '$Cumulative\: Number\:of\: Threads$')
plt.plot(dd_date, dd.year, 'r-o', label = '$Yearly\:Number\:of\:Threads$')
plt.legend(loc=2,numpoints=1,fontsize=13)
plt.show()

"""
Explanation: values_counts
End of explanation
"""


dg = data.groupby('year').sum()
dg

dgs = dg.cumsum()
dgs

def getDate(dat):
    dat_date_str = map(lambda x: str(x) +'-01-01', dat.index)
    dat_date = pd.to_datetime(dat_date_str)
    return dat_date

dg.date = getDate(dg)

fig = plt.figure(figsize=(12,5))
plt.plot(dg.date, dg.click, 'r-o', label = '$Yearly\:Number\:of\:Clicks$')
plt.plot(dg.date, dg.reply, 'g-s', label = '$Yearly\:Number\:of\:Replies$')
plt.plot(dg.date, dg.fans_num, 'b->', label = '$Yearly\:Number\:of\:Fans$')

plt.yscale('log')

plt.legend(loc=4,numpoints=1,fontsize=13)
plt.show()

data.groupby('year')['click'].sum()

data.groupby('year')['click'].mean()

"""
Explanation: groupby
End of explanation
"""


repost = []
for i in df.title:
    if u'转载' in i.decode('utf8'):
        repost.append(1)
    else:
        repost.append(0)

df['repost'] = repost

df.groupby('repost').sum()

df['click'][df['repost']==0][:5]

df['click'][df['repost']==1][:5]

from scipy import stats
stats.ttest_ind(data.click, data.repost)

sm.stats.ttest_ind(data.click, data.reply)
# test statistic, pvalue and degrees of freedom

"""
Explanation: 常用的统计分析方法

t检验
卡方检验
相关
回归

T-test
http://statsmodels.sourceforge.net/devel/stats.html
RQ: 转载的文章的点击量是否显著地小于原创的文章？
End of explanation
"""


from scipy.stats import chisquare
chisquare([16, 18, 16, 14, 12, 12], f_exp=[16, 16, 16, 16, 16, 8])


"""
Explanation: A chi-squared test
https://en.wikipedia.org/wiki/Chi-squared_test

also referred to as χ² test (or chi-square test), is any statistical hypothesis test in which the sampling distribution of the test statistic is a chi-square distribution when the null hypothesis is true. 
A chi-squared test can then be used to reject the null hypothesis that the data are independent.
Test statistics that follow a chi-squared distribution arise from an assumption of independent normally distributed data, which is valid in many cases due to the central limit theorem. 

Chi-squared tests are often constructed from a sum of squared errors, or through the sample variance.


Suppose there is a city of 1 million residents with four neighborhoods: A, B, C, and D. 


A random sample of 650 residents of the city is taken and their occupation is recorded as "blue collar", "white collar", or "no collar". 


The null hypothesis is that each person's neighborhood of residence is independent of the person's occupational classification. The data are tabulated as:


|               |  A           |    B          |C              | D  | Total|
| -------------|:-------------:|:-------------:|:-------------:|-----:|-----:|
| White collar| 90 | 60 | 104 |95 | 349|
| Blue collar| 30 | 50 | 51 | 20| 151|
| No coloar| 30 | 40 | 45 | 35|150|
| Total | 150 | 150| 200| 150| 650|

Let us take the sample living in neighborhood A, 150/650, to estimate what proportion of the whole 1 million people live in neighborhood A. 
Similarly we take 349/650 to estimate what proportion of the 1 million people are white-collar workers. 
By the assumption of independence under the hypothesis we should "expect" the number of white-collar workers in neighborhood A to be

$
\frac{150}{650} \frac{349}{650} 650 = 80.54
$
Then in that "cell" of the table, we have
$\frac{(\text{observed}-\text{expected})^2}{\text{expected}} = \frac{(90-80.54)^2}{80.54}$.
The sum of these quantities over all of the cells is the test statistic. 
Under the null hypothesis, it has approximately a chi-square distribution whose number of degrees of freedom are
$ (\text{number of rows}-1)(\text{number of columns}-1) = (3-1)(4-1) = 6. $
If the test statistic is improbably large according to that chi-square distribution, then one rejects the null hypothesis of independence.

scipy.stats.chisquare(f_obs, f_exp=None, ddof=0, axis=0)[source]

Calculates a one-way chi square test.
The chi square test tests the null hypothesis that the categorical data has the given frequencies.

Parameters: 
- f_obs : array_like Observed frequencies in each category.
- f_exp : array_like, optional Expected frequencies in each category. By default the categories are assumed to be equally likely.
- ddof : int, optional
End of explanation
"""


from scipy.stats import chisqprob, chi2
# p_value = chi2.sf(chi_statistic, df)
print chisqprob(3.94,1), 1 - chi2.cdf(3.94,1)


"""
Explanation: 
End of explanation
"""


print np.corrcoef(data.click, data.reply)

data.corr()


plt.plot(df.click, df.reply, 'r-o')
plt.show()

plt.plot(df.click, df.reply, 'gs')
plt.xlabel('$Clicks$', fontsize = 20)
plt.ylabel('$Replies$', fontsize = 20)
plt.xscale('log')
plt.yscale('log')
plt.title('$Allowmetric\,Law$', fontsize = 20)
plt.show()

"""
Explanation: Correlation
End of explanation
"""


import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Load data
dat = sm.datasets.get_rdataset("Guerry", "HistData").data
# Fit regression model (using the natural log of one of the regressors)
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()

"""
Explanation: Regression
End of explanation
"""


# Inspect the results
print results.summary()

reg = smf.ols('reply ~ click + followed_num', data=data).fit()

print reg.summary()

reg1 = smf.ols('np.log(reply+1) ~ np.log(click+1) +np.log(followed_num+1)+month', data=data).fit()
print reg1.summary()

fig = plt.figure(figsize=(12,8))
fig = sm.graphics.plot_partregress_grid(reg1, fig = fig)
plt.show()


import statsmodels.api as sm
from statsmodels.formula.api import ols

moore = sm.datasets.get_rdataset("Moore", "car",
                                 cache=True) # load data
data = moore.data
data = data.rename(columns={"partner.status" :
                             "partner_status"}) # make name pythonic
 



data[:5]

moore_lm = ols('conformity ~ C(fcategory, Sum)*C(partner_status, Sum)',
                 data=data).fit()

fig = plt.figure(figsize=(12,8))
fig = sm.graphics.plot_partregress_grid(moore_lm, fig = fig)
plt.show()

table = sm.stats.anova_lm(moore_lm, typ=2) # Type 2 ANOVA DataFrame
print table

"""
Explanation: 有些使用windows的同学无法运行上述代码：
在spyder中打开terminal
输入: pip install -U patsy
https://groups.google.com/forum/#!topic/pystatsmodels/KcSzNqDxv-Q
End of explanation
"""


