text:
Table of Contents

1  NEEDS TO BE UPDATED IN ACCORDANCE WITH ETF Market Internals_Test Class Script.ipynb

## NEEDS TO BE UPDATED IN ACCORDANCE WITH `ETF Market Internals_Test Class Script.ipynb`

```python
"""
Script to examine performance of Long/Short Portfolio

Factor: Cross-Sectional Return Bins

Strategy Classification: Long/Short Market-Neutral Momentum

1) Calculate return bins for each ETF Symbol's components (weekly)
2) Compare return bins across ETFs
3) Long N top ETFs that exceed PERCENTAGE threshold; Short N worst ETFs that exceed PERCENTAGE threshold.
4) Equal Weight, Risk-Parity, Minimum Correlation Algorithm
5) Hold N weeks and repeat
"""
```

    "\nScript to examine performance of Long/Short Portfolio\n\nFactor: Cross-Sectional Return Bins\n\nStrategy Classification: Long/Short Market-Neutral Momentum \n\n1) Calculate return bins for each ETF Symbol's components (weekly)\n2) Compare return bins across ETFs\n3) Long N top ETFs that exceed PERCENTAGE threshold; Short N worst ETFs that exceed PERCENTAGE threshold.\n4) Equal Weight, Risk-Parity, Minimum Correlation Algorithm\n5) Hold N weeks and repeat\n"

```python
%load_ext watermark
%watermark

%load_ext autoreload
%autoreload 2
# import standard libs
from IPython.display import display
from IPython.core.debugger import set_trace as bp
from pathlib import PurePath, Path
import sys
import time
from collections import OrderedDict as od
import re
import os
import json
os.environ['THEANO_FLAGS'] = 'device=cpu,floatX=float32'

# get project dir
pp = PurePath(Path.cwd()).parts[:-1]
pdir = PurePath(*pp)
script_dir = pdir / 'src'
viz_dir = pdir / 'report' / 'figures'
data_dir = pdir / 'data'
script_dirs = list(Path(script_dir).glob('./*/'))
for sdir in script_dirs: sys.path.append(sdir.as_posix())

import _blk_utilities_
import _blk_ETF_Internals_
utils = _blk_utilities_.utilities()

from decimal import Decimal
dplaces = Decimal('0.0001')

import pandas as pd
from pandas.tseries.offsets import *
import pandas_datareader.data as web
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
sz = 11
mlt = 1.66
size=(mlt * sz, sz)
import seaborn as sns
sns.set_style('white', {"xtick.major.size": 2, "ytick.major.size": 2})
flatui = ["#9b59b6", "#3498db", "#95a5a6", "#e74c3c", "#34495e", "#2ecc71","#f4cae4"]
sns.set_palette(sns.color_palette(flatui,7))
from collections import OrderedDict as od
from tqdm import tqdm, tqdm_notebook
import warnings
warnings.filterwarnings("ignore")
import missingno as msno

print()
%watermark -p pandas,pandas_datareader,dask,numpy,matplotlib,seaborn
```

    Wall time: 1.1 s

```python
def _quantize(val):
    return Decimal(val).quantize(dplaces)

def _format_bins(npos, nneg):
    npos['Percent of Total Stocks'] = npos['Percent of Total Stocks'].apply(_quantize)
    nneg['Percent of Total Stocks'] = nneg['Percent of Total Stocks'].apply(_quantize)
    return npos, nneg

"""Create combined dataframe"""
def _add_columns(df, empty_df, etf_sym):
    empty_df[etf_sym] = df['Percent of Total Stocks']
    return empty_df
```

```python
"""
import sample holdings from ETF
"""
def _remove_cash_holdings(etf):
    return pd.Series([sym for sym in etf if 'CASH' not in sym])
# ============================================================================================ #
# EUROPE
def _format_FTSE_symbols():
    ftse100 = pd.read_csv(project_path + r'_ETF_Holdings\_Foreign\ISF_holdings.csv', skiprows=2)
    ftse100 = ftse100[ ftse100['Sector'] != r'Cash and/or Derivatives' ].dropna()
    ftse100 = (ftse100['Issuer Ticker'].replace(['/', '\.'], '', regex=True)).sort_values().reset_index(drop=True)
    idx_loc = None
    for i in range(len(ftse100)):
        if ftse100.ix[i] == 'BTA':
            idx_loc = i
    ftse100[idx_loc] = 'BT-A'
    ftse100 =   pd.Series( [sec + '.L' for sec in ftse100 ] ).map(str.strip)
    return ftse100

def _format_DAX_symbols():
    dax = pd.read_csv(project_path + r'_ETF_Holdings\_Foreign\EXS1_holdings.csv', skiprows=2)
    dax = dax[ dax['Sector'] != r'Cash and/or Derivatives' ].dropna()
    dax = pd.Series( [sec + '.DE' for sec in dax['Issuer Ticker'] ] ).map(str.strip)
    return dax
# ============================================================================================ #
# CANADA

def _format_TSX_symbols():
    tsx = pd.read_csv(project_path + r'_ETF_Holdings\_Foreign\XIU_holdings.csv', skiprows=10)
    tsx = tsx[ tsx['Sector'] != r'Cash and/or Derivatives' ].dropna()
    tsx = utils._clean_symbols_strings(tsx['Ticker'])
    tsx = pd.Series( [sec + '.TO' for sec in tsx] ).map(str.strip)
    return tsx
# ============================================================================================ #
# AUSTRALIA
def _format_ASX_symbols():
    asx = pd.read_excel(project_path + r'_ETF_Holdings\_Foreign\fund_allholdings_STW.xls', skiprows=4)
    asx = asx[asx['Sector Classification'] != 'Unassigned']
    asx = asx['Ticker'].map(str.strip)
    asx = asx.replace('-AU', '.AX', regex=True)
    return asx
# ============================================================================================ #
# US MAJOR AVERAGES

def _format_SPY_symbols():
    spy = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-spy.xls', skiprows=3)
    spy = spy[spy['Shares Held']>0]
    spy = spy['Identifier'].map(str.strip)
    spy = utils._clean_symbols_strings(spy)
    spy = _remove_cash_holdings(spy)
    return spy

def _format_DIA_symbols():
    dia = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-dia.xls', skiprows=3)
    dia = dia[(dia['Shares Held']>0) & (dia['Identifier'] != 'CASH_USD')]
    dia = dia['Identifier'].map(str.strip)
    return dia

def _format_QQQ_symbols():
    qqq = pd.read_csv(project_path + r'_ETF_Holdings\_QQQ_holdings\holdings-QQQ.csv')
    qqq = qqq['HoldingsTicker'].map(str.strip)
    qqq = _remove_cash_holdings(qqq)
    return qqq
# ============================================================================================ #
# SECTOR SPDRS

def _format_XLB_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlb.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLE_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xle.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLF_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlf.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLI_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xli.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLK_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlk.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLP_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlp.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLU_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlu.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLV_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xlv.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

def _format_XLY_symbols():
    syms = pd.read_excel(project_path + r'_ETF_Holdings\_SPDR_holdings\holdings-xly.xls', skiprows=1)
    syms = syms['Symbol'].map(str.strip)
    syms = utils._clean_symbols_strings(syms)
    syms = _remove_cash_holdings(syms)
    return syms

# ============================================================================================ #
'''
def _format_STOXX_symbols():
    stoxx50 = pd.read_csv(project_path + r'_ETF_Holdings\_Foreign\EUE_holdings.csv', skiprows=2)
    stoxx50= stoxx50[ stoxx50['Sector'] != r'Cash and/or Derivatives' ].dropna()
    stoxx50 = pd.Series( [sec + '.PA' for sec in stoxx50['Issuer Ticker'] ] ).map(str.strip)
    return stoxx50

'''

# ============================================================================================ #
# EUROPE
FTSE = _format_FTSE_symbols()
DAX = _format_DAX_symbols()
# ============================================================================================ #
# CANADA
TSX = _format_TSX_symbols()
# ============================================================================================ #
# AUSTRALIA
ASX = _format_ASX_symbols()
# ============================================================================================ #
# US MAJOR AVERAGES

SPY = _format_SPY_symbols()
DIA = _format_DIA_symbols()
QQQ = _format_QQQ_symbols()
# ============================================================================================ #
# SECTOR SPDRS

XLB = _format_XLB_symbols()
XLE = _format_XLE_symbols()
XLF = _format_XLF_symbols()
XLI = _format_XLI_symbols()
XLK = _format_XLK_symbols()
XLP = _format_XLP_symbols()
XLU = _format_XLU_symbols()
XLV = _format_XLV_symbols()
XLY = _format_XLY_symbols()

all_etf_symbols = od([
                    ('SPY',SPY), ('DIA',DIA), ('QQQ',QQQ),
                    ('XLB',XLB), ('XLE',XLE), ('XLF',XLF),
                    ('XLI',XLI), ('XLK',XLK), ('XLP',XLP),
                    ('XLU',XLU), ('XLV',XLV), ('XLY',XLY),
                    ('FTSE',FTSE), ('DAX',DAX), ('TSX',TSX), ('ASX',ASX)
                    ])

```

```python
import importlib
importlib.reload(_blk_ETF_Internals_)
etfi = _blk_ETF_Internals_.etfInternals()

# get query dates anchored to last 52 weeks Monday
start_query, end_query = etfi._get_query_dates()
date_format = utils._create_date_format(start_query, end_query)
print(start_query, end_query, date_format)
```

    'creating directory...'
    'directory created [complete]'
    2014-03-31 00:00:00 2016-03-25 00:00:00 [Monday, 03.31.2014 - Friday, 03.25.2016]

```python
%%time
# Get ETF Indices price data for performance analytics
etf_syms = ['SPY', 'DIA', 'QQQ', 'XLB', 'XLE', 'XLF', 'XLI', 'XLK', \
            'XLP', 'XLU', 'XLV', 'XLY', '^FTSE', '^GDAXI', '^GSPTSE', '^AXJO']

def _get_px(symbol, start, end):
    return web.DataReader(symbol, 'yahoo', start, end)['Adj Close']

ETF_px = pd.DataFrame()
for sym in etf_syms:
    ETF_px[sym] = _get_px(sym, start_query, end_query)

rets = np.log(ETF_px / ETF_px.shift(1)).dropna()
print(rets.describe())
```

                  SPY         DIA         QQQ         XLB         XLE         XLF  \
    count  447.000000  447.000000  447.000000  447.000000  447.000000  447.000000
    mean     0.000413    0.000366    0.000671    0.000165   -0.000424    0.000198
    std      0.009223    0.009032    0.010998    0.011641    0.015942    0.010795
    min     -0.043019   -0.037366   -0.044690   -0.042811   -0.066400   -0.043664
    25%     -0.004124   -0.004256   -0.004531   -0.006514   -0.008980   -0.005332
    50%      0.000585    0.000743    0.000768    0.000822   -0.000258    0.000834
    75%      0.005132    0.005392    0.006418    0.006836    0.008051    0.006334
    max      0.037675    0.038670    0.049135    0.035693    0.048277    0.041351

                  XLI         XLK         XLP         XLU         XLV         XLY  \
    count  447.000000  447.000000  447.000000  447.000000  447.000000  447.000000
    mean     0.000380    0.000666    0.000641    0.000723    0.000431    0.000564
    std      0.009993    0.010557    0.007828    0.009859    0.010916    0.009985
    min     -0.036209   -0.040232   -0.035413   -0.042058   -0.043609   -0.040760
    25%     -0.005053   -0.004669   -0.004043   -0.004953   -0.005439   -0.004697
    50%      0.000541    0.000932    0.000898    0.001438    0.000970    0.001108
    75%      0.005837    0.006394    0.004968    0.006765    0.006722    0.006480
    max      0.030759    0.049169    0.029113    0.026591    0.041363    0.037271

                ^FTSE      ^GDAXI     ^GSPTSE       ^AXJO
    count  447.000000  447.000000  447.000000  447.000000
    mean    -0.000265    0.000036   -0.000134   -0.000077
    std      0.010541    0.013927    0.008903    0.009912
    min     -0.047798   -0.048165   -0.031745   -0.041765
    25%     -0.005741   -0.008695   -0.005202   -0.005716
    50%      0.000110    0.000746    0.000394    0.000440
    75%      0.005385    0.007712    0.004699    0.005897
    max      0.034971    0.048521    0.028963    0.029913
    Wall time: 3.71 s

```python
%%time
"""
Get all ETF Internals via Loop
"""
# ======================================================================= #
def _single_run_download(ETF_symbol, symbols, *args):
    """Get Symbol data"""
    location_h5 = etfi._get_symbol_px_data(ETF_symbol, symbols)
# ======================================================================= #
def _calc_return_thresholds(list_of_thresholds, log_returns, N_stocks, pos=None):
        """Function to calculate the quantity of stocks meeting or exceeding a return threshold

        Parameters:
        ===========
        list_of_thresholds: a list of each return level
        log_returns: dataframe of equity returns
        pos = boolean, if True calculates positive thresholds else negative

        Returns:
        ========
        N_thresh = pd.Series() with stock counts for each return threshold
        """
        if pos:
            N_thresh = pd.DataFrame(columns=['Count', 'Percent of Total Stocks'], index=list_of_thresholds)
            for thresh in list_of_thresholds:
                l = [1 for i in log_returns.ix[-1] if i > thresh]
                N_thresh.loc[thresh]['Count'] = sum(l)
            N_thresh['Percent of Total Stocks'] = (N_thresh['Count'] / N_stocks)
            return N_thresh
        else:
            N_thresh = pd.DataFrame(columns=['Count', 'Percent of Total Stocks'], index=list_of_thresholds)
            for thresh in list_of_thresholds:
                l = [1 for i in log_returns.ix[-1] if i <= thresh]
                N_thresh.loc[thresh]['Count'] = sum(l)
            N_thresh['Percent of Total Stocks'] = (N_thresh['Count'] / N_stocks)
            return N_thresh

# ======================================================================= #
def _get_return_bins(ETF_symbol, symbols, location_h5, start_date, end_date):
    """Extract Close Prices"""
    px = etfi._extract_adjclose_px(location_h5, symbols)

    """Set key parameters"""
    N_stocks = len(px.columns)
    log_returns = etfi._calc_log_returns(px)
    log_returns = log_returns.ix[start_date:end_date]

    """Calculate thresholds"""
    pos_thresh = [0., .01, .02, .03, .04, .05, .06, .07, .1, .15, .2]
    neg_thresh = [ flt * -1. for flt in pos_thresh ]

    N_pos = _calc_return_thresholds(pos_thresh, log_returns.resample('1W', how='sum'), N_stocks, pos=True)
    N_neg = _calc_return_thresholds(neg_thresh, log_returns.resample('1W', how='sum'), N_stocks)
    return N_pos, N_neg
# ======================================================================= #
def _weekly_returns(ETF_symbol, returns, start_date, end_date):
    """Set Key parameters"""
    selected_returns = returns.ix[start_date:end_date]
    by_week_returns = selected_returns.sum()
    return by_week_returns
# ======================================================================= #
def _convert_symbol(symbol, reverse=None):
    '''function to convert foreign indices symbol to yahoo finance symbol'''
    if not reverse:
        incompatible_syms = ['FTSE','DAX','TSX','ASX']
        compatible_syms = ['^FTSE', '^GDAXI', '^GSPTSE', '^AXJO']

        compatible_symbol = None # Declare symbol
        # iterate to check if any symbols are incompatible
        for inc in incompatible_syms:
            if symbol == inc:
                idx_loc = incompatible_syms.index(inc)
                compatible_symbol = compatible_syms[idx_loc]
                return compatible_symbol
        # if compatible_symbol is None then symbol is good; return it
        if not compatible_symbol:
            return symbol
    elif reverse:
        compatible_syms = ['FTSE','DAX','TSX','ASX']
        incompatible_syms = ['^FTSE', '^GDAXI', '^GSPTSE', '^AXJO']

        compatible_symbol = None # Declare symbol
        # iterate to check if any symbols are incompatible
        for inc in incompatible_syms:
            if symbol == inc:
                idx_loc = incompatible_syms.index(inc)
                compatible_symbol = compatible_syms[idx_loc]
                return compatible_symbol

        # if compatible_symbol is None then symbol is good; return it
        if not compatible_symbol:
            return symbol
# ======================================================================= #
```

    Wall time: 0 ns

```python
today = pd.datetime.today().date()
date_str = today.strftime('%m-%d-%y')

start = today - 252 * BDay()
idx_monday = pd.bdate_range(start, today - 5 * BDay(), freq='W-MON')
idx_friday = pd.bdate_range(start, today, freq='W-FRI')
print(idx_friday[:5])
print(idx_monday[:5])

dates = []
n = 0
for ifri in idx_friday[1:]:
    try:
        n += 1
        imon = idx_monday[n - 1]
        dates.append((ifri, imon))
    except: break
```

    DatetimeIndex(['2015-04-10', '2015-04-17', '2015-04-24', '2015-05-01',
                   '2015-05-08'],
                  dtype='datetime64[ns]', freq='W-FRI')
    DatetimeIndex(['2015-04-13', '2015-04-20', '2015-04-27', '2015-05-04',
                   '2015-05-11'],
                  dtype='datetime64[ns]', freq='W-MON')

```python
%%time
# ~ 7-8 minutes to run
idx = range(len(dates))
by_date_pos = od()
by_date_neg = od()

for i in idx:
    pos_df = pd.DataFrame()
    neg_df = pd.DataFrame()
    for etf_sym, symbols in all_etf_symbols.items():
        today = pd.datetime.today().date()
        date_str = today.strftime('%m-%d-%y')
        location = price_path + r'{}\{} ETF Market Internals Price Data {}.h5'.format(date_str, etf_sym, date_str)

        MON, FRI = 1, 0
        start_date = dates[i][MON]
        end_date = dates[i][FRI]

        npos, nneg = _get_return_bins(etf_sym, symbols, location, start_date, end_date)
        npos, nneg = _format_bins(npos, nneg)

        pos_df = _add_columns(npos, pos_df, etf_sym)
        neg_df = _add_columns(nneg, neg_df, etf_sym)

    by_date_pos[start_date] = pos_df
    by_date_neg[start_date] = neg_df
```

    Wall time: 7min 28s

```python
%%time
IDX = range(len(dates))
MON, FRI = 1, 0
by_date_returns = od()

for n in IDX:
    for sym in all_etf_symbols.keys():
        today = pd.datetime.today().date()
        start_date = dates[n][1] # [MON]
        end_date = dates[n][0] # [FRI]

        SYM = _convert_symbol(sym)
        tmp_rets = _weekly_returns(SYM, rets, start_date, end_date)

        idx = []
        for i in tmp_rets.index:
            sec = _convert_symbol(i, reverse=True)
            idx.append(sec)
        tmp_df = pd.DataFrame(index=idx)
        tmp_rets.index = idx
        tmp_df['weekly_returns'] = tmp_rets
        # print(tmp_df)

    by_date_returns[start_date] = tmp_rets
```

    Wall time: 668 ms

```python
"""
Get time series of returns for Symbols in the top quintile of return distribution

1. Select positive(negative) return bins
2. Group by Percent Total of return distro
3. Select group
4. Examine next week's returns
5. Store results
"""

# POSITIVE THRESHOLDS
mondays = [date[MON] for date in dates]
fridays = [date[FRI] for date in dates]
#len(list(pos_bin_dict.items())) # <= index due based on positive or negative thresholds exceeded
# ======================================================================= #
pos_thresh = [0., .01, .02, .03, .04, .05, .06, .07, .1, .15, .2]
neg_thresh = [ flt * -1. for flt in pos_thresh ]

pbt_dict = {} # agg pos_bin_dict by return thresholds
# calculate quintiles
for thresh in pos_thresh:
    try:
        pos_bin_dict = {}
        for k, v in by_date_pos.items():
            tmp = pd.DataFrame()
            tmp['Percent_Total'] = v.loc[thresh].astype(float)
            tmp['Quintiles'] = pd.qcut(tmp['Percent_Total'], 2, labels=False)
            pos_bin_dict[k] = tmp.sort_values('Percent_Total', ascending=False)
    except: pass
    pbt_dict[thresh] = pos_bin_dict

# ======================================================================= #
HPR = {} # hpr by thresh
for thresh in pbt_dict.keys():
    hpr_long = pd.Series() # long portfolio holding period rets
    hpr_short = pd.Series()
    LS = pd.Series()
    for k, v in pbt_dict[thresh].items():
        tgt_date = pd.to_datetime(k)
        # select quintile and get symbols
        top_q_syms = list(v[v['Quintiles']==0].index)
        bottom_q_syms = list(v[v['Quintiles']==1].index)
        try:
            for i, n in enumerate(dates):
                if n[MON] == tgt_date: # get location of date
                    next_monday_close = dates[i + 1][MON]
        except: continue

        hpr_long.loc[tgt_date.date().strftime('%m-%d-%y')] = by_date_returns[next_monday_close][top_q_syms].mean()
        hpr_short.loc[tgt_date.date().strftime('%m-%d-%y')] = by_date_returns[next_monday_close][bottom_q_syms].mean()

        long = hpr_long.loc[tgt_date.date().strftime('%m-%d-%y')]
        short = hpr_short.loc[tgt_date.date().strftime('%m-%d-%y')]

        LS.loc[tgt_date.date().strftime('%m-%d-%y')] = long - short
    HPR[thresh] = LS

# ======================================================================= #
for thresh in sorted(list(pbt_dict.keys())):
    try:
        print('positive threshold: {:3.2%} | mean return: {:3.2%} | sum: {:3.2%}'.format(thresh,
                                                                          Decimal(HPR[thresh].mean()).quantize(dplaces),
                                                                          Decimal(HPR[thresh].sum()).quantize(dplaces)
                                                                                ))
    except: pass
```

    positive threshold: 0.00% | mean return: 0.15% | sum: 7.57%
    positive threshold: 1.00% | mean return: 0.39% | sum: 7.10%
    positive threshold: 2.00% | mean return: 0.40% | sum: 7.23%
    positive threshold: 3.00% | mean return: 0.34% | sum: 6.14%
    positive threshold: 4.00% | mean return: 0.25% | sum: 4.55%
    positive threshold: 5.00% | mean return: 0.23% | sum: 1.39%
    positive threshold: 6.00% | mean return: 0.14% | sum: 0.82%
    positive threshold: 7.00% | mean return: 0.55% | sum: 2.19%
    positive threshold: 10.00% | mean return: NaN% | sum: 0.00%
    positive threshold: 15.00% | mean return: NaN% | sum: 0.00%
    positive threshold: 20.00% | mean return: NaN% | sum: 0.00%

```python
# plot cumulative returns from distro
for thresh in sorted(list(pbt_dict.keys())):
    try:
        fig, ax = plt.subplots()
        HPR[thresh].cumsum().plot(figsize=size, ax=ax, color=flatui[4])
        plt.title('Cumulative Mean Returns at threshold = {}'.format(thresh))
        #plt.axhline(0, lw=1, color='k')
    except: pass
```

```python
# NEGATIVE THRESHOLDS
mondays = [date[MON] for date in dates]
fridays = [date[FRI] for date in dates]
#len(list(pos_bin_dict.items())) # <= index due based on positive or negative thresholds exceeded
# ======================================================================= #
pos_thresh = [0., .01, .02, .03, .04, .05, .06, .07, .1, .15, .2]
neg_thresh = [ flt * -1. for flt in pos_thresh ]

nbt_dict = {} # agg pos_bin_dict by return thresholds
# calculate quintiles
for thresh in neg_thresh:
    try:
        neg_bin_dict = {}
        for k, v in by_date_neg.items():
            tmp = pd.DataFrame()
            #print(thresh)
            tmp['Percent_Total'] = v.loc[thresh].astype(float)
            tmp['Quintiles'] = pd.qcut(tmp['Percent_Total'], 2, labels=False)
            neg_bin_dict[k] = tmp.sort_values('Percent_Total', ascending=False)
            #print(tmp.sort_values('Percent_Total', ascending=False))
    except: pass
    nbt_dict[thresh] = neg_bin_dict
# ======================================================================= #

HPR_n = {} # hpr by thresh
for thresh in nbt_dict.keys():
    hpr_long = pd.Series() # long portfolio holding period rets
    hpr_short = pd.Series()
    LS = pd.Series()
    for k, v in nbt_dict[thresh].items():
        tgt_date = pd.to_datetime(k)
        # select quintile and get symbols
        top_q_syms = list(v[v['Quintiles']==1].index)
        bottom_q_syms = list(v[v['Quintiles']==0].index)
        try:
            for i, n in enumerate(dates):
                if n[MON] == tgt_date: # get location of date
                    next_monday_close = dates[i + 1][MON]
        except: continue

        hpr_long.loc[tgt_date.date().strftime('%m-%d-%y')] = by_date_returns[next_monday_close][top_q_syms].mean()
        hpr_short.loc[tgt_date.date().strftime('%m-%d-%y')] = by_date_returns[next_monday_close][bottom_q_syms].mean()

        long = hpr_long.loc[tgt_date.date().strftime('%m-%d-%y')]
        short = hpr_short.loc[tgt_date.date().strftime('%m-%d-%y')]

        LS.loc[tgt_date.date().strftime('%m-%d-%y')] = long - short
    HPR_n[thresh] = LS

# ======================================================================= #
for thresh in sorted(list(nbt_dict.keys()), reverse=True):
    try:
        print('negative threshold: {:3.2%} | mean return: {:3.2%} | sum: {:3.2%}'.format(thresh,
                                                                          Decimal(HPR_n[thresh].mean()).quantize(dplaces),
                                                                          Decimal(HPR_n[thresh].sum()).quantize(dplaces)
                                                                                ))
    except: pass
```

    negative threshold: -0.00% | mean return: 0.18% | sum: 8.76%
    negative threshold: -1.00% | mean return: 0.17% | sum: 8.43%
    negative threshold: -2.00% | mean return: 0.21% | sum: 7.63%
    negative threshold: -3.00% | mean return: 0.16% | sum: 5.59%
    negative threshold: -4.00% | mean return: 0.05% | sum: 1.62%
    negative threshold: -5.00% | mean return: 0.10% | sum: 2.14%
    negative threshold: -6.00% | mean return: NaN% | sum: 0.00%
    negative threshold: -7.00% | mean return: NaN% | sum: 0.00%
    negative threshold: -10.00% | mean return: NaN% | sum: 0.00%
    negative threshold: -15.00% | mean return: NaN% | sum: 0.00%
    negative threshold: -20.00% | mean return: NaN% | sum: 0.00%

```python
# plot cumulative returns from distro
for thresh in sorted(list(nbt_dict.keys()), reverse=True):
    try:
        fig, ax = plt.subplots()
        HPR_n[thresh].cumsum().plot(figsize=size, ax=ax, color=flatui[4])
        plt.title('Cumulative Mean Returns at threshold = {}'.format(thresh))
        plt.axhline(0, lw=1, color='k')
    except: pass
```

```python
# create portfolios based on pct pos, pct neg
import sympy

```

max_stars_repo_path: notebooks/ETF Component Return Bins Backtesting.ipynb


text:
# Error propagation with Sympy

Sympy is a Python module for symbolic computation (like Mathematica and Matlab) with an elegant Python design.

It can easily do mechanical tasks such as computing first derivatives for error propagation. We use this here to write a generic error propagation function in only 20+ lines.

```python
import sympy
print(f"Sympy version {sympy.__version__}")
```

    Sympy version 1.8

## Independent variables

```python
def value_and_covariance_gen(expr, variables):
    expr = sympy.parse_expr(expr)

    symbols = sympy.symbols(variables)
    cov_symbols = sympy.symbols(tuple("C_" + k for k in variables))
    expr2 = sum(expr.diff(s) ** 2 * c for s, c in zip(symbols, cov_symbols))
    expr2 = expr2.simplify() # recommended for speed and accuracy

    fval = sympy.lambdify(symbols, expr)
    fcov = sympy.lambdify(symbols + cov_symbols, expr2)

    def fn(**kwargs):
        x = tuple(v[0] for v in kwargs.values())
        c = tuple(v[1] for v in kwargs.values())
        return fval(*x), fcov(*x, *c)

    return fn

def value_and_covariance(expr, **kwargs):
    return value_and_covariance_gen(expr, tuple(kwargs))(**kwargs)
```

That's all, folks!

`value_and_covariance_gen` generates a Python function that computes the value of the expression and the propagated covariance.

`value_and_covariance` is just a shortcut to generate and immediately call the function. If the generated function is called with different values several times it is more efficient to generate it once with the first command and then call it several times.

Generating the propagating function takes about 6ms on my computer. The generated Python code is as fast to evaluate as any other function that calls numpy functions.

**Limitations of this implemention**
- Only independent (uncorrelated) inputs are supported here
- Expression must consist of basic math (must be parsable by Sympy)

### Example usage

```python
value_and_covariance("a + b", a=(1, 0.1), b=(2, 0.2)
```

    (3, 0.30000000000000004)

```python
value_and_covariance("s / (s + b)", s=(5, 0.5), b=(10, 0.1))
```

    (0.3333333333333333, 0.001037037037037037)

ext: ipynb

max_stars_repo_name: HDembinski/essays

max_issues_repo_name: HDembinski/essays

max_forks_repo_name: HDembinski/essays


text:
# PageRank: Eigenvectores y Eigenvalores

## 1. Introducción

Las búsquedas por Internet cobran cada vez más y más relevancia, por lo que los motores de búsqueda han evolucionado enormemente para atender la siempre creciente cantidad de información disponible en la red.  Uno de los precursores de los motores de búsqueda por Internet, es el algoritmo PageRank de Google.

En este trabajo se aborda el algoritmo de PageRank en su versión más simple, para entender los conceptos que lo conforman y cómo es que asigna la importancia a cada una de las páginas del resultado de una búsqueda en internet.  A este algoritmo se le conoce coloquialmente como *The 25 Billion Dollar Eigenvector*.  Revisemos los conceptos que lo conforman.

## 2. Teoría

### 2.1 Breve introducción a Teoría de Grafos

Un **grafo** es un par ordenado de dos conjuntos: nodos y vínculos, que se denota $G=(N,V)$. Un vínculo ${xy}$ se dice que conecta a los nodos $x$ y $y$, y es denominado $xy$. El orden de un grafo es igual a la cantidad de nodos existentes.

Un **grafo orientado** se define como aquel en el que todos los vínculos que lo componen tienen una orientación definida. Teniendo dos nodos $x$ y $y$ cualesquiera, debe verificarse que debe existir como mínimo un vínculo dirigido entre ellos. Básicamente, son grafos en los cuales se ha añadido una orientación a las aristas, representada gráficamente por una flecha.

El **grafo de la web** describe los enlaces dirigidos entre páginas de la red de Internet. Es un grafo orientado cuyos nodos corresponden a las páginas de Internet, y un vínculo orientado conecta la página $X$ con la página $Y$ si existe un hipervínculo en la página $X$ que refiere a la página $Y$.

#### 2.1.1 Representación de grafos

Existen diferentes formas de representar un grafo (simple), además de la geométrica y muchos métodos para almacenarlos en una computadora. La estructura de datos usada depende de las características del grafo y el algoritmo usado para manipularlo. Entre las estructuras más sencillas y usadas se encuentran las listas y las matrices, aunque frecuentemente se usa una combinación de ambas. Las listas son preferidas en grafos dispersos porque tienen un eficiente uso de la memoria. Por otro lado, las matrices proveen acceso rápido, pero pueden consumir grandes cantidades de memoria. Por el objetivo de este trabajo se hablará sobre la estructura matricial.

**Estructura matricial**
* Matriz de adyacencia: El grafo está representado por una matriz cuadrada M de tamaño $n^{2}$, donde $n$ es el número de nodos. Si hay un vínculo entre un nodo $x$ y un nodo $y$, entonces el elemento $m_{x,y}$ es $1$, de lo contrario, es $0$.
* Matriz de incidencia: El grafo está representado por una matriz de $V$ (vínculos) por $N$ (nodos), donde [nodo, vínculos] contiene la información del vínculo ($1$ si es conectado, $0$ si es no conectado).

### 2.2 ¿Qué es el PageRank?

El algoritmo PageRank es utilizado por Google Search para ayudarle a determinar la importancia o relevancia de una página web basándose en el grafo de la web. Fue desarrollado por los fundadores de Google, Larry Page y Sergey Brin, en la Universidad de Stanford mientras estudiaban el posgrado en ciencias de la computación.

El algoritmo PageRank se basa en la idea de que se puede juzgar la importancia de una página web mirando las páginas que contienen un vínculo hacia la misma. Para esto, se definen $O_i$ al conjunto de páginas que son referenciadas por $i$ (llamadas `outlinks`) e $I_i$ como el conjunto de páginas que hacen referencia a la página $i$ (llamadas `inlinks`). Si una página $A$ contiene un *outlink* a otra página $B$ se interpreta que la página $A$ considera que el contenido de $B$ es relevante para la temática abordada en $A$. Si existen muchos *inlinks* hacia $B$ se considera que es de común acuerdo que la página $B$ es importante. Por otro lado, si la página $B$ tiene solamente un *inlink* pero éste proviene de una página $C$ con autoridad (importancia) decimos que $C$ transfiere su autoridad a $B$, es decir, indica que $B$ es importante. Utilizando estos conceptos de importancia y autoridad el algoritmo PageRank asigna un rango a cada página basándose en sus *inlinks*.

No es el único algoritmo utilizado para el motor de búsqueda de Google, pero fue el primero y el más conocido.

### 2.3 Algoritmo

El algoritmo inicial del PageRank lo podemos encontrar en el documento original donde sus creadores presentaron el prototipo de Google: “The Anatomy of a Large-Scale Hypertextual Web Search Engine":

$$
PR(A) = (1-d) + d\sum_{i=1}^n \frac{PR(i)}{C(i)}
$$

Donde:
- $PR(A)$ es el PageRank de la página A.
- $d$ es un factor de amortiguación que tiene un valor entre 0 y 1.
- $PR(i)$ son los valores de PageRank que tienen cada una de las páginas $i$ que enlazan a $A$.
- $C(i)$ es el número total de enlaces salientes de la página $i$ (sean o no hacia $A$).

#### 2.3.1 Factor de amortiguación

La teoría del PageRank sostiene que un usuario cualquiera que aleatoriamente hace click en distintos links eventualmente dejará de hacerlo. La probabilidad, en cualquier paso, de que la persona continúe es un factor de amortiguación $d$. Este modelo se puede entender como una cadena de Markov, que es un proceso estocástico donde el valor presente de una variable aleatoria es relevante para predecir el valor futuro de dicha variable. Esto nos indica que la probabilidad de realización de determinado suceso sólo depende de lo ocurrido en el periodo inmediatamente anterior.

Varios estudios han probado diferentes factores de amortiguación, pero generalmente se asume que el factor de amortiguación se establecerá alrededor de $0.85$.

La introducción del factor de amortiguación en la fórmula resta algo de peso a todas las páginas de Internet y consigue que las páginas que no tienen *outlinks* no salgan especialmente beneficiadas. Si un usuario aterriza en una página sin *outlinks*, lo que hará será navegar a cualquier otra página aleatoriamente, lo que equivale a suponer que una página sin *outlinks* tiene enlaces a todas las páginas de Internet.

### 2.4 Fórmula iterativa

$$
PR_{t+1}(P_i) = \sum_{P_j}\frac{PR_t(P_j)}{C(P_j)}
$$

Donde:
* $PR_{t+1}(P_i)$: PageRank de la página $i$ en la siguiente iteración.
* $PR_{t}(P_j)$: PageRank de la página $j$ en la iteración previa.
* $C(P_j)$: número total de *outlinks* de la página $j$.

En la iteración $0$, se inicializa un vector de tamaño $n$, y sus valores inciales para cada posición serán $\frac{1}{n}$.

### 2.5 Representación en matriz

Se pueden utilizar operaciones matriciales en vez de realizar la forma iterativa explicada previamente (donde se hacen cálculos y actualizaciones de valores de uno por uno). Con las operaciones matriciales se pueden realizar múltiples cálculos al mismo tiempo.

Se define la matriz $H$ que será la matriz de adyacencia asociada al grafo. Esta matriz se definirá por los *outlinks* y el número de éstos por cada página. La matriz $H$ será **estocástica** que nos ayudará a describir las transiciones en una cadena de Markov.

El método utilizado para realizar estos cálculos de forma más sencilla es el **método de la potencia**.

#### 2.5.1 Método de la potencia

$$ PR_{t+1} = H PR_t$$

Donde:
* $H$ es la matriz estocástica.

Esta matriz se debe multiplicar con un vector en cada iteración, ¿pero cuál es el vector inicial? Como se ha mencionado previamente, el vector inicial será un vector de tamaño $n$, y sus valores inciales para cada posición serán $\frac{1}{n}$.

Entonces, el vector de la siguiente iteración esta dado por:
    $$v_2 = H v$$

Y la siguiente iteración será:

$$v_3 = H v_2 = H (H v) = H^2 v$$

Entonces:

$$v_n = H^{(n-1)} v$$

Si hacemos varias iteraciones, se converge a un valor de equilibrio, el cual será el eigenvector asociado al eigenvalor con módulo máximo.

Esta es la principal diferencia del método de la potencia respecto a los algoritmos QR y de rotaciones de Jacobi, que aproximan a todos los eigenvalores y eigenvectores.  Adicionalmente, al sólo requerirse una operación matriz-vector, resulta más eficiente para matrices grandes y ralas, como las resultantes de una búsqueda de internet (que en el 2006, tenían una dimensión de aproximadamente $8.1\times10^9$ x $8.1\times10^9$).

## 3. Datos

Para fines ilustrativos y mostrar con claridad el funcionamiento del algoritmo, se utilizará un set de datos sencillo en el reporte.  Si se desea ahondar con sets de datos que contengan más nodos y vínculos, se puede hacer referencia al [Ejemplo de implementación](https://github.com/lecepe00/PageRank-eigenvectores-eigenvalores/tree/main/Implementacion) incluido en el repositorio.

## 4. Implementación

### Ejemplo:

Supongamos que tenemos 4 páginas web donde:
- la página A apunta a las páginas C y B
- la página B apunta a la página D
- la página C apunta a las páginas A, B, D
- la página D apunta a la página C

Veámoslo de forma gráfica:

```python
import numpy as np
import math
import pandas as pd
from scipy.sparse import csr_matrix

# Para las pruebas
from pytest import approx
import pprint

# Para los grafos
import networkx as nx
import matplotlib.pyplot as plt
```

```python
pages_links = pd.DataFrame( columns=['pagina', 'outlinks'])
pages_links['pagina'] = ['A', 'B', 'C', 'D']
pages_links['outlinks'][0] = ['C', 'B']
pages_links['outlinks'][1] = ['D']
pages_links['outlinks'][2] = ['A', 'B', 'D']
pages_links['outlinks'][3] = ['C']
```

```python
nodos = pages_links['pagina']

G = nx.DiGraph()

for nodo in nodos:
    G.add_node(nodo)

for page, outlinks in zip(pages_links['pagina'],pages_links['outlinks']):
    for outlink in outlinks:
        G.add_edge(page, outlink)

pos = nx.spring_layout(G)

nx.draw_networkx_nodes(G, pos, node_color='#a0a88a', node_size=400)
nx.draw_networkx_labels(G, pos)
nx.draw_networkx_edges(G, pos, edge_color='#305029', width=1.8, arrows = True, style='dashed')

plt.show()
```

Definimos nuestra matriz estocástica conforme a lo previamente definido:

```python
sparseMatrix = csr_matrix((len(pages_links), len(pages_links)),
                          dtype = np.float).toarray()
newdf = pd.DataFrame(sparseMatrix, columns = pages_links['pagina'], index=pages_links['pagina'])
count = 0
for page, outlinks in zip(pages_links['pagina'],pages_links['outlinks']):
    score = 1/len(outlinks)
    dict_page = {page: np.zeros(len(pages_links))}
    for outlink in outlinks:
        newdf[page][outlink] = score
```

```python
H = newdf.to_numpy()
H
```

    array([[0.        , 0.        , 0.33333333, 0.        ],
           [0.5       , 0.        , 0.33333333, 0.        ],
           [0.5       , 0.        , 0.        , 1.        ],
           [0.        , 1.        , 0.33333333, 0.        ]])

Ahora definimos el vector inicial, que en este caso será:

\begin{align}
    v &= \begin{bmatrix}
           1/4 \\
           1/4 \\
           1/4 \\
           1/4
         \end{bmatrix}
  \end{align}

```python
v = np.ones(len(pages_links))
v = v*1/len(v)
pprint.pprint(v)
```

    array([0.25, 0.25, 0.25, 0.25])

Como habíamos definido antes:

$$v_n = H^{n-1}v$$

donde $n$ será el número de iteraciones.

En este caso probaremos con 20 iteraciones para obtener el vector 21.

```python
n = 20
H = np.linalg.matrix_power(H, n)
v_n = H@v
```

Cuando hacemos los cálculos obtenemos $v_n$:

```python
v_n
```

    array([0.12500368, 0.18750089, 0.37499297, 0.31250246])

Para mejor entendimiento:

```python
df = pd.DataFrame(v_n, columns=['pagerank'])

df['nodo'] = pages_links['pagina']
df = df.set_index('nodo')

df
```

      pagerank

      nodo

      A
      0.125004

      B
      0.187501

      C
      0.374993

      D
      0.312502

Volvemos a graficar con sus respectivos PageRank y podemos observar que $C$ es la página con más alto PageRank, por lo tanto será mostrado primero en el motor de búsqueda.

```python
power_method = df.pagerank.to_dict()

power_method_round = {}
for node in power_method:
    if node in power_method:
        #set the node name as the key and the label as its value
        power_method_round[node] = round(power_method[node],2)
```

```python
G = nx.DiGraph()

for nodo in power_method:
    G.add_node(nodo)

for page, outlinks in zip(pages_links['pagina'],pages_links['outlinks']):
    for outlink in outlinks:
        G.add_edge(page, outlink)

pos = nx.spring_layout(G)

labels = {}
for node in G.nodes():
    if node in power_method:
        #definimos el nombre del nodo como la llave y la etiqueta como el valor
        labels[node] = "{}:{}".format(node,round(power_method[node],2))

# Es lo mismo que power_method pero redondeado
print(labels)

nx.draw_networkx_nodes(G, pos, node_color='#a0a88a', node_size=400)
nx.draw_networkx_labels(G, pos, labels=labels, font_weight="bold")
nx.draw_networkx_edges(G, pos, edge_color='#305029', width=1.8, arrows = True, style='dashed')

plt.show()
```

## 5. Conclusiones

El algoritmo de PageRank, aunque incorpora conceptos relativamente sencillos de grafos, álgebra lineal y probabilidad, los integra de una manera muy inteligente, que ayuda a darle orden a la enorme cantidad de información disponible en la red.  Es esta capacidad y flexibilidad de entender una base de datos descentralizada, siempre cambiante y en continuo crecimiento, lo que ha mantenido vigente a este algoritmo, a más de 20 años de su publicación original.

## 6. Referencias

- American Mathematical Society.  [How Google Finds Your Needle in the Web's Haystack](http://www.ams.org/publicoutreach/feature-column/fcarc-pagerank)
- Barriola, Juan M. y Dotta, Milena.  2015.  [¿CÓMO FUNCIONA GOOGLE? EL ALGORITMO PAGERANK, DIAGRAMAS DE GRAFOS Y CADENAS DE MARKOV](http://www.economicas.uba.ar/wp-content/uploads/2016/04/2-%C2%BFC%C3%B3mo-funciona-google_-el-algoritmo-pagerank-diagramas-de-grafos-y-cadenas-de-markov.-Juan-Manuel-Barriola-y-Milena-Dotta-2.pdf)
- Brin, Sergey & Page, Lawrence.  [The Anatomy of a Large-Scale Hypertextual Web Search Engine](https://web.archive.org/web/20170624024436/http://infolab.stanford.edu/~backrub/google.html)
- Bryan, Kurt & Leise, Tanya.  [THE LINEAR ALGEBRA BEHIND GOOGLE](https://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf)
- Equipo 4.  Optimización 2020-2.  [Ejemplo de implementación](https://github.com/lecepe00/PageRank-eigenvectores-eigenvalores/tree/main/Implementacion)
- GeeksforGeeks.  [PageRank](https://www.geeksforgeeks.org/page-rank-algorithm-implementation/)
- Global Software Support.  [PageRank Algorithm - Power Method](https://video.search.yahoo.com/yhs/search?fr=yhs-ddc-linuxmint&hsimp=yhs-linuxmint&hspart=ddc&p=page+rank#id=9&vid=b35a1c4b4379bde2aeb371cd71f68b6b&action=view)
- Ken State University.  [Eigenvalue Computations](http://www.math.kent.edu/~reichel/courses/intr.num.comp.2/lecture21/evmeth.pdf)
- Palacios M., Erick. (2020).  Libro de Optimización 2020.  [Aplicaciones y usos: Page Ranking en un buscador de páginas web](https://itam-ds.github.io/analisis-numerico-computo-cientifico/II.computo_matricial/2.3/Algoritmos_y_aplicaciones_de_eigenvalores_eigenvectores_de_una_matriz.html#aplicaciones-y-usos-page-ranking-en-un-buscador-de-paginas-web)
- Wikipedia.  [PageRank](https://en.wikipedia.org/wiki/PageRank)
- Wikipedia.  [Teoría de grafos](https://es.wikipedia.org/wiki/Teor%C3%ADa_de_grafos)

```python

```


text:
## Maximum Likelihood Estimation

Maximum likelihood estimation is one of the key techniques employed in statistical signal processing for a wide variety of applications from signal detection to parameter estimation. In the following, we consider a simple experiment and work through the details of maximum likelihood estimation to ensure that we understand the concept in one of its simplest applications.

### Setting up the Coin Flipping Experiment

Suppose we have coin and want to estimate the probability of heads ($p$) for it. The coin is Bernoulli distributed:

$$ \phi(x)= p^x (1-p)^{(1-x)} $$

where $x$ is the outcome, *1* for heads and *0* for tails. The $n$ independent flips, we have the likelihood:

$$ \mathcal{L}(p|\mathbf{x})= \prod_{i=1}^n p^{ x_i }(1-p)^{1-x_i} $$

This is basically notation. We have just substituted everything into $ \phi(x)$ under the independent-trials assumption.

The idea of *maximum likelihood* is to maximize this as the function of $p$ after plugging in all of the $x_i$ data. This means that our estimator, $\hat{p}$ , is a function of the observed $x_i$ data, and as such, is a random variable with its own distribution.

### Simulating the Experiment

We need the following code to simulate coin flipping.

```
from __future__ import division
from scipy.stats import bernoulli
import numpy as np

p_true=1/2 # this is the value we will try to estimate from the observed data
fp=bernoulli(p_true)

def sample(n=10):
    'simulate coin flipping'
    return fp.rvs(n)# flip it n times

xs = sample(100) # generate some samples
```

Now, we can write out the likelihood function using `sympy`

```
import sympy
from sympy.abc import x, z
p=sympy.symbols('p',positive=True)

L=p**x*(1-p)**(1-x)
J=np.prod([L.subs(x,i) for i in xs]) # objective function to maximize
```

Below, we find the maximum using basic calculus. Note that taking the `log` of $J$ makes the maximization problem tractable but doesn't change the extrema.

```
logJ=sympy.expand_log(sympy.log(J))
sol=sympy.solve(sympy.diff(logJ,p),p)[0]

x=linspace(0,1,100)
plot(x,map(sympy.lambdify(p,logJ,'numpy'),x),sol,logJ.subs(p,sol),'o',
                                          p_true,logJ.subs(p,p_true),'s',)
xlabel('$p$',fontsize=18)
ylabel('Likelihood',fontsize=18)
title('Estimate not equal to true value',fontsize=18)
```

Note that our estimator $\hat{p}$ (red circle) is not equal to the true value of $p$ (green square), but it is at the maximum of the likelihood function. This may sound disturbing, but keep in mind this estimate is a function of the random data; and since that data can change, the ultimate estimate can likewise change. I invite you to run this notebook a few times to observe  this. Remember that the estimator is a *function* of the data and is thus also a *random variable*, just like the data is.

Let's write some code to empirically examine the behavior of the maximum likelihood estimator using a simulation of multiple trials. All we're doing here is combining the last few blocks of code.

```
def estimator_gen(niter=10,ns=100):
    'generate data to estimate distribution of maximum likelihood estimator'
    out=[]
    x=sympy.symbols('x',real=True)
    L=   p**x*(1-p)**(1-x)
    for i in range(niter):
        xs = sample(ns) # generate some samples from the experiment
        J=np.prod([L.subs(x,i) for i in xs]) # objective function to maximize
        logJ=sympy.expand_log(sympy.log(J))
        sol=sympy.solve(sympy.diff(logJ,p),p)[0]
        out.append(float(sol.evalf()))
    return out if len(out)>1 else out[0] # return scalar if list contains only 1 term

etries = estimator_gen(100) # this may take awhile, depending on how much data you want to generate
hist(etries) # histogram of maximum likelihood estimator
title('$\mu=%3.3f,\sigma=%3.3f$'%(mean(etries),std(etries)),fontsize=18)
```

Note that the mean of the estimator ($\mu$) is pretty close to the true value, but looks can be deceiving. The only way to know for sure is to check if the estimator is unbiased, namely, if

$$ \mathbb{E}(\hat{p}) = p $$

Because this problem is simple, we can solve for this in general noting that since $x=0$ or $x=1$, the terms in the product of $\mathcal{L}$ above are either $p$, if $x_i=1$ or $1-p$ if $x_i=0$. This means that we can write

$$ \mathcal{L}(p|\mathbf{x})= p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i} $$

with corresponding log as

$$ J=\log(\mathcal{L}(p|\mathbf{x})) =  \log(p)  \sum_{i=1}^n x_i +   \log(1-p) \left(n-\sum_{i=1}^n x_i\right)$$

Taking the derivative of this gives:

$$  \frac{dJ}{dp} = \frac{1}{p}\sum_{i=1}^n x_i + \frac{(n-\sum_{i=1}^n x_i)}{p-1} $$

and solving this leads to

$$  \hat{p} = \frac{1}{ n} \sum_{i=1}^n x_i $$

This is our *estimator* for $p$. Up til now, we have been using `sympy` to solve for this based on the data $x_i$ but now we have it generally and don't have to solve for it again. To check if this estimator is biased, we compute its expectation:

$$ \mathbb{E}\left(\hat{p}\right) =\frac{1}{n}\sum_i^n \mathbb{E}(x_i) = \frac{1}{n} n \mathbb{E}(x_i) $$

by linearity of the expectation and where

$$\mathbb{E}(x_i)  = p$$

Therefore,

$$ \mathbb{E}\left(\hat{p}\right) =p $$

This means that the esimator is unbiased. This is good news. We almost always want our estimators to be unbiased. Similarly,

$$ \mathbb{E}\left(\hat{p}^2\right) = \frac{1}{n^2} \mathbb{E}\left[\left(  \sum_{i=1}^n x_i \right)^2 \right]$$

and where

$$ \mathbb{E}\left(x_i^2\right) =p$$

and by the independence assumption,

$$ \mathbb{E}\left(x_i x_j\right) =\mathbb{E}(x_i)\mathbb{E}( x_j) =p^2$$

Thus,

$$ \mathbb{E}\left(\hat{p}^2\right) =\left(\frac{1}{n^2}\right) n
\left[
p+(n-1)p^2
\right]
$$

So, the variance of the estimator, $\hat{p}$ is the following:

$$ \sigma_\hat{p}^2 = \mathbb{E}\left(\hat{p}^2\right)- \mathbb{E}\left(\hat{p}\right)^2  = \frac{p(1-p)}{n} $$

Note that the $n$ in the denominator means that the variance asymptotically goes to zero as $n$ increases (i.e. we consider more and more samples). This is good news also because it means that more and more coin flips leads to a better estimate of the underlying $p$.

Unfortunately, this formula for the variance is practically useless because we have to know $p$ to compute it and $p$ is the parameter we are trying to estimate in the first place! But, looking at $ \sigma_\hat{p}^2 $, we can immediately notice that if $p=0$, then there is no estimator variance because the outcomes are guaranteed to be tails. Also, the maximum of this variance, for whatever $n$, happens at $p=1/2$. This is our worst case scenario and the only way to compensate is with more samples (i.e. larger $n$).

All we have computed is the mean and variance of the estimator. In general, this is insufficient to characterize the underlying probability density of $\hat{p}$, except if we somehow knew that  $\hat{p}$ were normally distributed. This is where the powerful [*central limit theorem*](http://mathworld.wolfram.com/CentralLimitTheorem.html) comes in. The form of the estimator, which is just a mean estimator, implies that we can apply this theorem and conclude that  $\hat{p}$ is normally distributed. However, there's a wrinkle here: the theorem tells us that  $\hat{p}$ is asymptotically normal, it doesn't quantify how many samples $n$ we need to approach this asymptotic paradise. In our simulation this is no problem since we can generate as much data as we like, but in the real world, with a costly experiment, each sample may be precious. In the following, we won't apply this theorem and instead proceed analytically.

### Probability Density for the Estimator

To write out the full density for $\hat{p}$, we first have to ask what is the probability that the estimator will equal a specific value and the tally up all the ways that could happen with their corresponding probabilities. For example, what is the probability that

$$  \hat{p} = \frac{1}{n}\sum_{i=1}^n x_i  = 0 $$

This can only happen one way: when $x_i=0 \hspace{0.5em} \forall i$. The probability of this happening can be computed from the density

$$ f(\mathbf{x},p)= \prod_{i=1}^n \left(p^{x_i} (1-p)^{1-x_i}  \right) $$

$$ f\left(\sum_{i=1}^n x_i  = 0,p\right)= \left(1-p\right)^n  $$

Likewise, if $\lbrace x_i \rbrace$ has one $i^{th}$  value equal to one, then

$$ f\left(\sum_{i=1}^n x_i  = 1,p\right)= n p \prod_{i=1}^{n-1} \left(1-p\right)$$

where the $n$ comes from the $n$ ways to pick one value equal to one from the $n$ elements $x_i$. Continuing this way, we can construct the entire density as

$$ f\left(\sum_{i=1}^n x_i  = k,p\right)= \binom{n}{k} p^k  (1-p)^{n-k}  $$

where the term on the left is the binomial coefficient of $n$ things taken $k$ at a time. This is the binomial distribution and it's not the density for $\hat{p}$, but rather for $n\hat{p}$. We'll leave this as-is because it's easier to work with below. We just have to remember to keep track of the $n$ factor.

#### Confidence Intervals

Now that we have the full density for $\hat{p}$, we are ready to ask some meaningful questions. For example,

$$ \mathbb{P}\left( | \hat{p}-p | \le \epsilon p \right) $$

Or, in words, what is the probability we can get within $\epsilon$ percent of the true value of $p$. Rewriting,

$$ \mathbb{P}\left(  p - \epsilon p \lt \hat{p} \lt p + \epsilon p \right) = \mathbb{P}\left(  n p - n \epsilon p \lt \sum_{i=1}^n x_i   \lt n p + n \epsilon p \right)$$

Let's plug in some live numbers here for our worst case scenario where $p=1/2$. Then, if $\epsilon = 1/100$, we have

$$ \mathbb{P}\left( \frac{99 n}{100} \lt \sum_{i=1}^n x_i   \lt \frac{101 n}{100} \right)$$

Since the sum in integer-valued, we need $n> 100$ to even compute this. Thus, if $n=101$ we have

$$ \mathbb{P}\left( \frac{9999}{200} \lt \sum_{i=1}^{101} x_i   \lt \frac{10201}{200} \right) = f\left(\sum_{i=1}^{101} x_i  = 50,p\right)= \binom{101}{50} (1/2)^{50}  (1-1/2)^{101-50} = 0.079$$

This means that in the worst-case scenario for $p=1/2$, given $n=101$ trials, we will only get within 1% of the actual $p=1/2$ about 8% of the time. If you feel disappointed, that only means you've been paying attention. What if the coin was really heavy and it was costly to repeat this 101 times? Then, we would be within 1% of the actual value only 8% of the time. Those odds are terrible.

Let's come at this another way: given I could only flip the coin 100 times, how close could I come to the true underlying value with high probability (say, 95%)? In this case we are seeking to solve for $\epsilon$. Plugging in gives,

$$ \mathbb{P}\left(  50 - 50 \epsilon \lt \sum_{i=1}^{100} x_i   \lt 50 + 50 \epsilon  \right) = 0.95$$

which we have to solve for $\epsilon$. Fortunately, all the tools we need to solve for this are already in `scipy`.

```
import scipy.stats

b=scipy.stats.binom(100,.5) # n=100, p = 0.5, distribution of the estimator \hat{p}

f,ax= subplots()
ax.stem(arange(0,101),b.pmf(arange(0,101))) # heres the density of the sum of x_i

g = lambda i:b.pmf(arange(-i,i)+50).sum() # symmetric sum the probability around the mean
print 'this is pretty close to 0.95:%r'%g(10)
ax.vlines( [50+10,50-10],0 ,ax.get_ylim()[1] ,color='r',lw=3.)

```

The two vertical lines in the plot show how far out from the mean we have to go to accumulate 95% of the probability. Now, we can solve this as

$$ 50 + 50 \epsilon = 60 $$

which makes $\epsilon=1/5$ or 20%. So, flipping 100 times means I can only get within 20% of the real $p$ 95% of the time in the worst case scenario (i.e. $p=1/2$).

```
b=scipy.stats.bernoulli(.5) # coin distribution
xs = b.rvs(100) # flip it 100 times
phat = mean(xs) # estimated p

print abs(phat-0.5) < 0.5*0.20 # did I make it w/in interval 95% of the time?
```

    True

Let's keep doing this and see if we can get within this interval 95% of the time.

```
out=[]
b=scipy.stats.bernoulli(.5) # coin distribution
for i in range(500): # number of tries
    xs = b.rvs(100) # flip it 100 times
    phat = mean(xs) # estimated p
    out.append(abs(phat-0.5) < 0.5*0.20 ) # within 20%

print 'Percentage of tries within 20 interval = %3.2f'%(100*sum(out)/float(len(out) ))
```

    Percentage of tries within 20 interval = 96.20

Well, that seems to work. Now we have a way to get at the quality of the estimator, $\hat{p}$.

## Summary

In this section, we explored the concept of maximum likelihood estimation using a coin flipping experiment both analytically and numerically with the  scientific Python tool chain. There are two key points to remember. First, maximum likelihood estimation produces a function of the data that is itself a random variable, with its own statistics and distribution. Second, it's worth considering how to analytically derive the density function of the estimator rather than relying on canned packages to compute confidence intervals wherever possible. This is especially true when data is hard to come by and the approximations made in the central limit theorem are therefore harder to justify.

### References

This [IPython notebook](www.ipython.org) is available for [download](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Maximum_likelihood.ipynb). I urge you to experiment with the calculations for different parameters. As always, corrections and comments are welcome!

max_stars_repo_name: fperez/Python-for-Signal-Processing

max_issues_repo_name: fengyemuse/Python-for-Signal-Processing

max_forks_repo_name: fengyemuse/Python-for-Signal-Processing


text:
# 8. Higher order derivatives

we consider the second derivative problem

$$
	a(v, u) = \int_{\Omega} \nabla^2 u  \nabla^2 v + u v
    ~d\Omega, \quad \forall u,v \in \mathcal{V}_h
$$

```python
import numpy as np
from numpy import linspace, zeros, pi

from sympy.core.containers import Tuple
from sympy import symbols
from sympy import Symbol
from sympy import Lambda
from sympy import Function

from gelato.glt import glt_symbol
from gelato.calculus   import (Dot, Cross, Grad, Curl, Rot, Div, dx)
from gelato.calculus   import Constant
from gelato.fem.assembly import assemble_matrix
from gelato.fem.utils    import compile_kernel
from gelato.fem.utils    import compile_symbol

from spl.fem.splines import SplineSpace
from spl.fem.tensor  import TensorFemSpace
from spl.fem.vector  import VectorFemSpace
```

```python
x = Symbol('x')

u = Symbol('u')
v = Symbol('v')

a = Lambda((x,v,u), dx(dx(u)) * dx(dx(v)) + u*v)
```

```python
# ...  create a finite element space
p  = 3
ne = 64

grid = linspace(0., 1., ne+1)

V = SplineSpace(p, grid=grid, nderiv=2)
# ...
```

```python
kernel_b = compile_kernel('kernel_b', a, V, backend='fortran')
```

```python
M = assemble_matrix(V, kernel_b).tocsr()
M = M.todense()
```

```python
# import eigenvalue solver
from scipy.linalg import eig

w, v = eig(M)
wr = w.real
wr.sort()
```

```python
expr = glt_symbol(a, space=V)
```

```python
symbol_b = compile_symbol('symbol_b', a, V, backend='fortran')
```

```python
t1 = linspace(-pi,pi, ne+1)
x1 = linspace(0.,1., ne+1)
e = zeros(ne+1)
symbol_b(x1, t1, e)
e.sort()
```

```python
%matplotlib inline
import matplotlib.pyplot as plt
```

```python
plt.plot(e, "+b", label="glt symbol")
plt.plot(wr, "xr", label="eigenvalues")
plt.legend(loc=2);
```

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:600px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    h2 {
        font-family: 'Fenix', serif;
    }
    h3{
		font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
       }
	h4{
		font-family: 'Fenix', serif;
       }
    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }
    div.text_cell_render{
        font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 135%;
        font-size: 120%;
        width:600px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro";
			font-size: 90%;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h1 {
        font-weight: 200;
        font-size: 50pt;
		line-height: 100%;
        color:#054BCD;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #054BCD;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>

```python

```

max_stars_repo_name: ratnania/IGA-Python

max_issues_repo_name: GabrielJie/IGA-Python

max_forks_repo_name: GabrielJie/IGA-Python


text:
# Introduction to Linear Algebra

## What is Linear Algebra?
Linear algebra is the branch of mathematics concerning linear equations such as linear maps such as and their representations in vector spaces and through matrices. Linear algebra is central to almost all areas of mathematics [See More](https://en.wikipedia.org/wiki/Linear_algebra)

## Applications of Linear Algebra in Data Science
- Loss Functions
- Regularization
- Covariance Matrix
- Support Vector Machine Classification
- Principal Component Analysis (PCA)
- Singular Value Decomposition
- Word Embeddings
- Latent Semantic Analysis (LSA)
- Image Representation as Tensors
- Convolution and Image Processing

__SEE DETAILS:__
https://www.analyticsvidhya.com/blog/2019/07/10-applications-linear-algebra-data-science/

## Linear Algebra Operations

```python
# import numpy
import numpy as np
```

### Determinant of matrix
- The determinant of a matrix is a special number that can be calculated from a **square matrix**
\begin{equation}
A = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \\
det(A) = ad - bc
\end{equation}
Whre, A is a $2 \times 2$ matrix.

- `np.linalg.det()` - performs determinant of the matrix
- **Syntax:** `np.linalg.det(matrix)`

__Source:__
https://www.mathsisfun.com/algebra/matrix-determinant.html

```python
A = np.arange(0, 16).reshape(4,4)
print(A)
```

    [[ 0  1  2  3]
     [ 4  5  6  7]
     [ 8  9 10 11]
     [12 13 14 15]]

```python
# determinant of A
np.linalg.det(A)
```

    -2.9582283945787796e-30

```python
B = np.matrix("4,5,6,7;2,-3,3,3; 3,4,5,6; 4, 7,8,9")
```

```python
# determinant of B
np.linalg.det(B)
```

    12.000000000000005

### Rank of matrix
- `np.linalg.matrix_rank()` - returns rank of the matrix
- **Syntax:** `np.linalg.matrix_rank(matrix)`

```python
# rank
np.linalg.matrix_rank(A)
```

    2

```python

```

```python

```

```python

```

```python

```

```python
### Rank of matrix
- `np.linalg.matrix_rank()` - returns rank of the matrix
- **Syntax:** `np.linalg.matrix_rank(matrix)`
```

```python

```

## References
- https://numpy.org/
- https://www.edureka.co/blog/python-numpy-tutorial/
- https://github.com/enthought/Numpy-Tutorial-SciPyConf-2019
- [Python Machine Learning Cookbook](https://www.amazon.com/Python-Machine-Learning-Cookbook-Prateek/dp/1786464470)
<hr>

*This notebook was created by [Jubayer Hossain](https://jhossain.me/) | Copyright © 2020, [Jubayer Hossain](https://jhossain.me/)*

max_stars_repo_name: hossainlab/dsnotes

max_issues_repo_name: hossainlab/dsnotes

max_forks_repo_name: hossainlab/dsnotes


text:
# Dimentionality Reduction

## PCA

### Maximum variance formulation

consider a dataset $\{x_{i}\}$ where $i=1,...,n$ and $x_{i} \in \mathbb{R}^{d}$.

our goal is to project the data onto a space having dimensionality $k < d$ while maximizing the variance of the projected data.

to begin with, consider the projection onto a one-dimensional space$(k=1)$.

we can define the direction of this space by a vector $u_{1} \in \mathbb{R}^{d}$, we can choose $u_{1}$ to be a unit vector so that $u_{1}^{T}u_{1} = 1$.

each data point $x_{i}$ is then projected onto a scalar value $u_{1}^{T}x_{i}$, then mean of the projected data:

$$\frac{1}{n}\sum_{i=1}^{n}u_{1}^{T}x_{i} = u_{1}^{T}\overline{x}$$

the variance of the projected data:

$$\frac{1}{n}\sum_{i=1}^{n}(u_{1}^{T}x_{i} - u_{1}^{T}\overline{x})^{2} = \frac{1}{n}\sum_{i=1}^{n}u_{1}^{T}(x_{i} - \overline{x})(x_{i} - \overline{x})^{T}u_{1} = u_{1}^{T}Su_{1}$$

where

$$S = \frac{1}{n}\sum_{i=1}^{n}(x_{i} - \overline{x})(x_{i} - \overline{x})^{T}$$

now we can formalize our problem as:

$$\underset{u_{1}}{min}\ -u_{1}^{T}Su_{1}$$
$$s.t\quad u_{1}^{T}u_{1} = 1$$

the lagrangian of this optimization problem:

$$L(u_{1}, \lambda_{1}) = -u_{1}^{T}Su_{1} + \lambda_{1}(u_{1}^{T}u_{1} - 1)$$

the primal:

$$\underset{u_{1}}{min}\ \underset{\lambda_{1}}{max}\ L(u_{1}, \lambda_{1})$$

primal satisfy the KKT conditions, so equivalent to dual:

$$\underset{\lambda_{1}}{max}\ \underset{u_{1}}{min}\ L(u_{1}, \lambda_{1})$$

setting the derivative with respect to $u_{1}$ equal to zero, we have:

$$Su_{1} = \lambda_{1}{u_{1}}$$

which say that $u_{1}$ must be a eigenvector of $S$, if we left-multiply by $u_{1}^{T}$ and make use of $u_{1}^{T}u_{1} = 1$, we get:

$$u_{1}^{T}Su_{1} = \lambda_{1}$$

and so the variance will be a maximum when we set $u_{1}$ equal to the eigenvector having the largest eigenvalue $\lambda_{1}$. this eigenvector is known as the first principal component.

we can define the additional principal components in an increamental fashion by choosing each new direction to be that which maximizes the projected variance amongst all possible directions orthogonal to those already considered.

second principal component:

$$\underset{u_{2}}{min}\ -u_{2}^{T}Su_{2}$$
$$s.t\quad u_{2}^{T}u_{2} = 1, u_{1}^{T}u_{2} = 0$$

like before, using lagrangian we derive:

$$Su_{2} = \lambda_{2}{u_{2}} + \phi{u_{1}}$$

left multiply by $u_{1}^{T}$:

$$u_{1}^{T}Su_{2} = \lambda_{2}u_{1}^{T}{u_{2}} + \phi{u_{1}^{T}}{u_{1}}$$

analyzing each component:

$$u_{1}^{T}Su_{2} = u_{2}^{T}Su_{1} = u_{2}^{T}\lambda_{1}u_{1} = \lambda{u_{1}^{T}u_{2}} = 0$$
$$u_{1}^{T}{u_{2}} = 0$$
$${u_{1}^{T}}{u_{1}} = 1$$

we get:

$$\phi = 0$$

back to zero derivative we have:

$$Su_{2} = \lambda_{2}{u_{2}}$$
$$u_{2}^{T}Su_{2} = \lambda_{2}$$

so $\lambda_{2}$ is the second largest eigenvalue of $S$.

by induction, we can show that $i$-th principal component is the $i$-th largest eigenvector of $S$.

### properties of non-negative definite symmetric real matrix

$$S = \frac{1}{n}\sum_{i=1}^{n}(x_{i} - \overline{x})(x_{i} - \overline{x})^{T}$$

is of that kind.

### Minimum-error formulation

a complete orthonormal basis vectors $u_{i}$ in $\mathbb{R}^{d}$:

$$u_{i}^{T}u_{j} = \delta_{ij}$$

$x_{k}$ coordinate with respect to $u_{i}$ is $x_{k}^{T}u_{i}$, so:

$$x_{k} = \sum_{i=1}^{d}(x_{k}^{T}u_{i})u_{i}$$

$x_{k}$ can be approximated by the $m$-dimensional subspace representation plus a constant:

$$\tilde{x}_{k} = \sum_{i=1}^{m}z_{ki}u_{i} + \sum_{i=m+1}^{d}b_{i}u_{i}$$

where $z_{ki}$ depend on the particular data point, whereas ${b_{i}}$ are constants that are the same for all data points.

our goal is to minimize:

$$J = \frac{1}{n}\sum_{k=1}^{d}\left \| x_{k} - \tilde{x}_{k} \right \|^{2} $$

setting the derivative with respect to $z_{ni}$ to zero, and making use of the orthonormality conditions, we obtain:

$$z_{ni} = x_{n}^{T}u_{i}$$

similarly, we obtain:

$$b_{i} = \overline{x}^{T}u_{i}$$

substitude for $z_{ni}$ and $b_{i}$, we obtain:

$$x_{k} - \tilde{x}_{k} = \sum_{i=m+1}^{d}((x_{k} - \overline{x}_{k})^{T}u_{i})u_{i}$$

finally our goal is to minimize:

$$J = \frac{1}{n}\sum_{k=1}^{n}\sum_{i=m+1}^{d}(x_{k} - \overline{x}_{k})^{2} = \sum_{i=m+1}^{d}u_{i}^{T}Su_{i}$$

this is similar to the maximum variance formulation in the opposite direction.

### manual data

```python
"""construct dataset"""
import numpy as np

m = 100
w1, w2 = 0.1, 0.3
noise = 0.2

X = np.empty((m, 3))
X[:, :2] = np.random.multivariate_normal([0, 0], [[2, 1], [1, 5]], m)
X[:, 2] = X[:, 0] * w1 + X[:, 1] * w2 + noise * np.random.randn(m)
```

```python
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
X2D = pca.fit_transform(X)

X2D[: 5]
```

    array([[-0.40565221, -1.69570572],
           [-1.65068551,  0.30621659],
           [-0.4562352 ,  0.33505227],
           [ 1.23654634, -1.16955408],
           [-2.39337295,  0.04917042]])

```python
"""sklearn actually uses SVD"""
X_centered = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(X_centered)

X2D_SVD = X_centered.dot(Vt.T[:, :2])
X2D_SVD[: 5]
```

    array([[ 0.40565221,  1.69570572],
           [ 1.65068551, -0.30621659],
           [ 0.4562352 , -0.33505227],
           [-1.23654634,  1.16955408],
           [ 2.39337295, -0.04917042]])

### mnist data

```python
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split

mnist = fetch_openml('mnist_784', version=1, as_frame=False)
mnist.target = mnist.target.astype(np.uint8)

X = mnist["data"]
y = mnist["target"]

X_train, X_test, y_train, y_test = train_test_split(X, y)
```

```python
"""
set n_components (0.0, 1.0), indicating the ratio of variance you wish to preserve
use inverse_transform trying to inverse
"""
pca = PCA(n_components=0.95)
X_reduced = pca.fit_transform(X_train)
X_mnist = pca.inverse_transform(X_reduced)
```

```python
pca.explained_variance_ratio_[: 10]
```

    array([0.09756334, 0.07156108, 0.06168814, 0.05388291, 0.04877275,
           0.04311341, 0.03261168, 0.0289652 , 0.02763173, 0.02347036])

```python
"""IncrementalPCA"""
from sklearn.decomposition import IncrementalPCA

n_batches = 100
inc_pca = IncrementalPCA(n_components=154)
for X_batch in np.array_split(X_train, n_batches):
    print(".", end="")
    inc_pca.partial_fit(X_batch)

X_reduced = inc_pca.transform(X_train)
```

    ....................................................................................................

### swiss roll

```python
from sklearn.datasets import make_swiss_roll

X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=42)
```

```python
"""kernel pca"""
from sklearn.decomposition import KernelPCA

rbf_pca = KernelPCA(n_components=2, kernel="rbf", gamma=0.04)
X_reduced = rbf_pca.fit_transform(X)
```

## Locally Linear Embedding(LLE)

LLE works by first measuring how each training instance linearly relates to it's colsest neighbors (c.n)

then looking for a low-dimensional representation of the training set where these local relationships are best preserved.

LLE step one: linearly modeling local relationships:

$$
\begin{equation}
\begin{split}
\hat{W} &= \underset{W}{argmin}\sum_{i=1}^{m}\left(x^{(i)} - \sum_{j=1}^{m}w_{ij}x^{(j)}\right)^{2}\\
\mbox{s.t }\ &1.w_{ij} = 0 \mbox{ if } x^{(j)} \mbox{ is not one of c.n of } x^{(i)} \\
&2.\sum_{i=1}^{m}w_{ij} = 1 \mbox{ for all }j
\end{split}
\end{equation}
$$

LLE second step doing the reverse: keeping the weights fixed and finding the optimal position of the instances' image in low-dimensional space, suppose $x^{(i)}$'s low-dimensional image is $z^{(i)}$.

$$
\begin{equation}
\begin{split}
\hat{Z} &= \underset{Z}{argmin}\sum_{i=1}^{m}\left(z^{(i)} - \sum_{j=1}^{m}w_{ij}z^{(j)}\right)^{2}\\
\mbox{s.t }\ &1.\sum_{i=1}^{m}z^{(i)} = 0\\
&2.\sum_{i=1}^{m}(z^{(i)}z^{(i)})^{T} = mI_{d}
\end{split}
\end{equation}
$$

```python
from sklearn.manifold import LocallyLinearEmbedding

lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10)
X_reduced = lle.fit_transform(X)
```

## Exercise

```python
from sklearn.ensemble import RandomForestClassifier
from sklearn.base import clone

rnd_clf = RandomForestClassifier(n_estimators=1000, max_depth=4, n_jobs=-1, random_state=42)
low_clf = clone(rnd_clf)
```

```python
from datetime import datetime
from sklearn.metrics import accuracy_score

pre_time = datetime.now()
rnd_clf.fit(X_train, y_train)
print((datetime.now() - pre_time).seconds)
accuracy_score(y_test, rnd_clf.predict(X_test))
```

    19

    0.8194285714285714

```python
pca = PCA(n_components=0.95)
X_train_reduced = pca.fit_transform(X_train)
X_test_reduced = pca.transform(X_test)
```

```python
pre_time = datetime.now()
low_clf.fit(X_train_reduced, y_train)
print((datetime.now() - pre_time).seconds)
accuracy_score(y_test, low_clf.predict(X_test_reduced))
```

    20

    0.8050857142857143

```python
np.random.seed(42)

m = 10000
idx = np.random.permutation(60000)[:m]

X = mnist['data'][idx]
y = mnist['target'][idx]
```

```python
from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, random_state=42)
X_reduced = tsne.fit_transform(X)
```

```python
import matplotlib.pyplot as plt

plt.figure(figsize=(13,10))
plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=y, cmap="jet")
plt.axis('off')
plt.colorbar()
plt.show()
```

```python

```

max_stars_repo_name: newfacade/machine-learning-notes

max_issues_repo_name: newfacade/machine-learning-notes

max_forks_repo_name: newfacade/machine-learning-notes